Gabriel Peyré
email: gabriel.peyre@ceremade.dauphine.fr

Best Basis Compressed Sensing

This paper proposes an extension of compressed sensing that allows to express the sparsity prior in a dictionary of bases. This enables the use of the random sampling strategy of compressed sensing together with an adaptive recovery process that adapts the basis to the structure of the sensed signal. A fast greedy scheme is used during reconstruction to estimate the best basis using an iterative refinement. Numerical experiments on sounds and geometrical images show that adaptivity is indeed crucial to capture the structures of complex natural signals.

Introduction

Classical Sampling vs. Compressed Sensing

The classical sampling theory of Shannon is based on uniform smoothness assumptions (low frequency spectral content). Under band limited condition, finely enough sampled functions can be recovered from a set of n pointwise measurements.

However, most natural signals f are characterized by very different prior assumptions, such as a decomposition with few elements on some fixed orthogonal basis B. This is the case for a sparse expansion of a sound in a 1D local Fourier basis or the compression of a natural image using a wavelet expansion. Under such sparseness assumption, one can hope to use a much smaller number n < N of measurements, which are linear projections Φf = { f, ϕ i } n i=1 on a set of fixed vectors ϕ i ∈ R N . The price to pay for this compressed sampling strategy is a non-linear reconstruction procedure to recover f from the compressed representation Φf . This theory of compressed acquisition of data has been pushed forward during last few years conjointly by Candès and Tao [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] and Donoho [START_REF] Donoho | Compressed sensing[END_REF].

In order for this recovery to be effective, one needs sensing vectors ϕ i that are incoherent with the vectors of B. A convenient way to achieve this property is to use random vectors ϕ i , which cannot be sparsely represented in basis B.

Application in imaging. Compressed sensing acquisition of data could have an important impact for the design of imaging devices where data acquisition is expensive. For instance in seismic or magnetic imaging one could hope to use few random projections of the object to acquire together with a high precision reconstruction.

Analogies in physiology. This compressed sampling strategy could potentially lead to interesting models for various sensing operations performed biologically. Skarda and Freeman [START_REF] Skarda | Does the brain make chaos in order to make sense of the world?[END_REF] have proposed a non-linear chaotic dynamic to explain the analysis of sensory inputs. This chaotic state of the brain ensures robustness toward unknown events and unreliable measurements, without using too many computing resources. While the theory of compressed sensing is presented here as a random acquisition process, its extension to deterministic or dynamic settings is a fascinating area for future research in signal processing.

The Best Basis Approach

Frames vs. dictionary of bases. Fixed orthogonal bases are not flexible enough to capture the complex redundancy of sounds or natural images. For instance the orthogonal wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] lacks of translation and rotation invariance and is not efficient to compress geometric images [START_REF] Donoho | Wedgelets: Nearly minimax estimation of edges[END_REF][START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF]. It is thus useful to consider families of vectors that are redundant but offer a stable decomposition. For instance, frames of translation invariant wavelets have been used for image denoising and frames of rotation-invariant Gabor functions are useful to characterize textures [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

However, to capture the complex structure of sounds or the geometry of natural images, one needs a very large number of such elementary atoms. Frame theory suffers from both theoretical difficulties (lack of stability) and technical problems (computational complexity) when the number of basis vectors increases too much. To cope with these problems, one can consider a dictionary D = {B λ } λ∈Λ of orthogonal bases B λ . Choosing an optimal basis in such a dictionary allows to adapt the approximation to the complex content of a specific signal.

Sound processing. Atoms with a broad range of localizations in space and frequency are needed to represent the transient features that exist in sounds such as the one depicted in figure 3, (a). Local cosine bases [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] divide the time axis in segments that are adapted to the local frequency content of the sound. Other kinds of dictionaries of 1D bases have been proposed, such as the wavelet packets dictionary [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF].

Images and geometry. The set of cartoon images is a simple model that captures the sketch content of natural images. Figure 4,(a), shows such a geometrically regular image, which contains smooth areas surrounded by regular curves. The curvelet frame of Candès and Donoho [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] can deal with such a regularity and enjoys a better approximation rate than traditional isotropic wavelets. This result can be enhanced using a dictionary of locally elongated functions that follow the image geometry. Bandelets bases of Le Pennec and Mallat [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF] provide such a geometric dictionary together with a fast optimization procedure to compute a basis adapted to a given image.

Adaptive biological computation. Hubel and Wiesel have shown that low level computation done in area V1 of the visual cortex are well approximated by multiscale oriented linear projections [START_REF] Hubel | Receptive fields and functional architecture of monkey striate cortex[END_REF]. Olshausen and Field proposed in [START_REF] Olshausen | Emergence of simple-cell receptive-field properties by learning a sparse code for natural images[END_REF] that redundancy is important to account for sparse representation of natural inputs. However further non-linear processings are done by the cortex to remove high order geometrical correlations present in natural images. Such computations are thought to perform long range groupings over the first layer of linear responses [START_REF] Lee | Computations in the early visual cortex[END_REF] and thus correspond to an adaptive modification of the overall neuronal response. The best-basis coding strategy could thus offer a signal-processing counterpart to this neuronal adaptivity. This paper only deals with orthogonal best basis search but extensions of this approach to dictionaries of redundant transforms are likely to improve numerical results and better cope with biological computations.

Best Basis Computation

Notations. The ℓ p norms are defined by

||f || p ℓ p = i |f [i]| p and ||f || ℓ 0 = # {i \ f [i] = 0} . A dictionary D Λ = {B λ } λ∈Λ is a set of orthogonal bases B λ = {ψ λ m } m of R N . A cost pen(λ)
is defined as a prior complexity measure associated to each basis B λ and satisfies λ∈Λ 2 -pen(λ) = 1. A fixed weight pen(λ) = log 2 (M ) can be used if the size M of D is finite. The parameter pen(λ) can be interpreted as the number of bits needed to specify a basis B λ . Following the construction of Coifman and Wickerhauser [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF], a best basis B λ * adapted to a signal minimizes a Lagrangian L

λ * = argmin λ∈Λ L(f, λ, t)
where

L(f, λ, t) = ||Ψ λ f || ℓ 1 + C 0 t pen(λ), (1)
where

Ψ λ = [ψ λ 0 , . . . , ψ λ N -1 ]
T is the transform matrix defined by B λ . This best basis B λ * is thus the one that gives the sparsest description of f as measured by the ℓ 1 norm. The Lagrange multiplier C 0 t weights the penalization pen(λ) associated to the complexity of a basis B λ . The parameter t is the noise level due to acquisition errors or approximate sparsity, whereas C 0 is a scaling constant that can be tuned for a specific dictionary D Λ . Note that t is estimated iteratively during the recovery algorithms presented in subsections 2.2 and 3.2.

Practical best basis computation. For a dictionary D that enjoys a multiscale structure, the optimization of L is carried out with a fast procedure. This algorithm uses the fact that if a basis splits as

B λ = B λ0 ∪ B λ1 , the Lagrangian satisfies L(f, λ, t) = L(f, λ 0 , t) + L(f, λ 1 , t).
It is thus only needed to compute the value of L on elementary bases from which each basis B λ can be decomposed. The number of such elementary bases is usually much smaller than the total number of bases. A bottom-up regression algorithm is then used to search for λ * , see [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF][START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF] for practical examples of this process.

Compressed Sensing Reconstruction

In this paper, the sampling matrix Φ = [ϕ 0 , . . . , ϕ n-1 ] T is defined by random points {ϕ i } i of unit length although other random sensing schemes can be used, see [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF].

Basis Pursuit Formulation

Searching for the sparsest signal f * in some basis B = {ψ m } m that matches the sensed values y = Φf leads to consider

f * = argmin g∈R N ||Ψ g|| ℓ 0 subject to Φg = y.
(

) 2
where

Ψ = [ψ 0 , . . . , ψ N -1 ]
T is the transform matrix defined by B. The combinatorial optimization ( 2) is NP-hard to solve and convexification of the objective function is introduced in the following basis pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] formulation

f * = argmin g∈R N ||Ψ g|| ℓ 1 subject to Φg = y. (3)
Candès and Tao in [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] and Donoho in [START_REF] Donoho | Compressed sensing[END_REF] have shown that, if f is sparse enough in some basis B, f can be recovered from the sensed data y = Φf . More precisely, they show that it exists a constant

C such that if ||Ψ f || ℓ 0 k then if n C log(N ) k one recovers f * = f .
To deal with noisy measurements y = Φf + w, where w is a Gaussian noise of variance t 2 , one can turn the constrained formulations (3) into a penalized variational problem

f * = argmin g∈R N 1 2 ||Φg -y|| 2 ℓ 2 + t||Ψ g|| ℓ 1 . (4)
The Lagrange multiplier t accounts both for stabilisation against noise and approximate sparsity, which is common in practical applications. This approximate sparsity is characterized by the decays of the best M -term approximation f M of f in basis B defined by

f M = m∈It f, ψ m ψ m where I t = {m \ | f, ψ m | > t}
and M = Card(I t ). This can also be understood as an approximation using a thresholding of the coefficients of f in B that are below t. The approximation power of B is measured by the decay of the approximation error ||f -f M || C M -α . The parameter α is a measure of the smoothness of f with respect to the basis B.

The function f with an approximation exponent α can be seen as a k-sparse function f corrupted by a deterministic noise of amplitude k -α . One thus needs to choose a threshold t ∼ p -α for a compressed sensing scenario where p C log(N ) k. Candès et al. [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF] shows that the ℓ 1 minimization (4) leads to a recovery error of order n -α up to some logarithmic factors.

The relationship t ∼ p -α (up to logarithmic factors) is not precise enough to be useful in real applications, but reveals the insight behind the regularization formulation [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. When the number n of acquired samples decreases, the noise in the function reconstructed with basis pursuit (3) increases and the regularization imposed by the sparsity ||Ψ f || ℓ 1 needs to be increased.

Iterative Thresholding For Sparsity Minimization

The recovery procedure suggested by equation ( 4) corresponds to the inversion of the operator Φ under sparsity constraints on the observed signal f . In this paper we follow Daubechies et al. that propose in [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] an algorithm based on iterated thresholding to perform such regularized inversion. This algorithm had been previously used for image restauration by Figueiredo and Nowak [START_REF] Figueiredo | An em algorithm for wavelet-based image restoration[END_REF] and is derived in [START_REF] Candès | Practical signal recovery from random projections[END_REF] as the iteration of two projections on convex sets.

Other greedy algorithms such as orthogonal matching pursuit (OMP) [START_REF] Tropp | Signal recovery from partial information via orthogonal matching pursuit[END_REF] and stagewise orthogonal matching pursuit [START_REF] Donoho | Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit[END_REF] (StOMP) have been used to solve the compressed sensing reconstruction in a fast way. These algorithms do not fit very well into the best-basis extension exposed in section 3, mainly because OMP is too computationally intensive for imaging applications and because StOMP performed a very agressive iterations (which speed up computation but bias the choice of best-basis toward the initially chosen basis).

We make repeated use of the soft thresholding operator defined in an orthogonal basis

B = {ψ m } m S t (g) = m s t ( g, ψ m ) ψ m where s t (x) = sign(x)(|x| -t) + .
The algorithm. The steps of the algorithm are • Initialization. Set s = 0, f 0 = 0.

• Step 1: Updating the estimate. Set f s+1 = f s + Φ T (y -Φ f s ). This update is an orthogonal projection on the convex set {g \ Φg = y}. • Step 2: Denoising the estimate. Estimate the noise level σ s using a median estimator in the transformed domain σ s = median(|Ψ f s+1 |)/0.6745 and set the threshold to t = 3σ s . Compute f s+1 = S t (f s+1 ) the thresholding of f in basis B. This step corresponds to a projection on a ℓ 1 ball {g \ ||Ψ g|| ℓ 1 c} for some c that depends on t. • Stopping criterion. If s < s max , go to step 1, otherwise stop iterations. In all our experiments, the number of iterations is set to s max = 20. The strategy of updating t through iterations has been adopted by Donoho et al. in their extension of orthogonal matching pursuit [START_REF] Donoho | Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit[END_REF]. The median estimator of the noise is based on the assumption, accurate in practice, that the current estimate f s+1 is corrupted by a Gaussian noise which diminishes during iterations.

Best Basis Compressed Sensing

Variational Formulation

To enhance the quality of the recovery in real life compressed sensing applications, the use of redundant frames has been proposed by several authors [START_REF] Donoho | Extensions of compressed sensing[END_REF][START_REF] Candès | Practical signal recovery from random projections[END_REF]. Redundancy can impact compressed sensing efficiency, because random sensing vectors can be correlated once inverted through the frame operator Ψ , since it is not orthogonal. This paper uses a different approach relying on orthogonal transforms which does not lower the recovery property of the sampling matrix Φ. It also allows a more efficient approximation through the use of a best basis in a large dictionary.

The compressed sensing machinery is extended to a dictionary of bases D Λ by imposing that the recovered signal is sparse in at least one basis of D Λ . To avoid using too complex basis the recovery process from noisy measurements takes into account a complexity pen(λ) of the optimal basis B λ . This is coherent with the best-basis approximation scheme of [START_REF] Donoho | Wedgelets: Nearly minimax estimation of edges[END_REF][START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF], although penalization is not strictly required in classical statistical estimation. The original recovery procedure ( 4) is replaced by the following minimization

f * = argmin g∈R N min λ∈Λ 1 2 ||Φg -y|| 2 ℓ 2 + t ||Ψ λ g|| ℓ 1 + C 0 t 2 pen(λ) , (5)
the penalization C 0 t 2 pen(λ) is the same as in equation (1).

Best Basis Recovery Algorithm

Searching in the whole dictionary D for the best basis that minimizes formulation ( 5) is not feasible for large dictionaries, which typically contain of the order of 2 N bases. Instead we propose a greedy search for the best basis during the recovery process. This leads to the following algorithm.

• Initialization. Set s = 0, f 0 = 0 and choose λ 0 ∈ Λ at random or using some default choice (such as a DCT basis in 1D or a wavelet basis in 2D).

• Step 1: Updating the estimate. Set f s+1 = f s + Φ T (y -Φ f s ).
• Step 2: Denoising the estimate. Compute the noise level

σ s = median(|Ψ f s+1 |)/0.6745
and set the threshold to t = 3σ s . Compute f s+1 = S t (f s+1 ) where S t is the threshold operator at t in the basis B λs . • Step 3: Update best basis. Compute λ s+1 = argmin λ L(f s+1 , λ, t). For typical dictionaries such as the ones considered in this paper, this minimization is carried out with a fast procedure, as seen in subsection 1.3. • Stopping criterion. If s < s max , go to step 1, otherwise stop iterations.

Adapted Local Cosine Transform

For each scale j > 0, the set of locations {0, . . . , N -1} is subdivided using N/2 j intervals [x j p , x j p+1 ], where the endpoints are given by x j p = 2 -j N p -1/2. A local cosine basis B p j = {ψ jp k } k is defined for each of these intervals using

∀ k = 0, . . . , N/2 j -1, ∀ ℓ, ψ jp k = b 2 j (ℓ -x j p ) 2 2 -j N cos π k + 1 2 ℓ -x j p 2 -j N ,
where b is a smooth windowing function that satisfies some compatibility conditions [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. A local cosine basis B λ of R N is parameterized by a binary tree segmentation λ of {0, . . . , N -1}, see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. The set of leaves L = {L p j } of λ are indexed by their depth j > 0 and position p in the binary tree, see figure 1, left. The basis B λ is the union of the various elementary bases B λ = B p j \ L p j ∈ L , see figure 1, right. For this local cosine basis dictionary, the penalization pen(λ) is defined as the number of leaves in the binary tree λ.

The decomposition of a signal f in each of the basis B p j is computed in O(N log(N ) 2 ) time using FFT. From these atomic decompositions, a best basis B λ * that minimizes (1) can be extracted by a tree pruning procedure in O(N ) time, see [START_REF] Coifman | Entropy-based algorithms for best basis selection[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

j =0 j = 1 j = 2 p = 1 p = 1 2 1 2 3 4
Fig. 1. A dyadic tree λ defining a spatial segmentation (left) ; some local cosine basis functions ψ jp k of the basis B λ (right).

Numerical Results

A synthetic sparse signal f = (Ψ λ ) -1 h is generated using a random local cosine basis B λ and a random signal of spikes h with ||h|| ℓ 0 = 30, see figure 2,(a). The signal recovered by the non-adaptative algorithm of subsection 2.2 in an uniform cosine basis B λ0 is significantly different from the original, figure 2,(b). This is due to the fact that f is less sparse in B λ0 , since

||Ψ λ0 f || ℓ 0 = 512 (a) (c1) (c2) (c3) (b)

and ||Ψ

λ0 f || ℓ 1 ≈ 2.8||Ψ λ f || ℓ 1 .
During the iterations of the algorithm presented in subsection 3.2, the estimated best basis B λs evolves in order to match the best basis B λ , see figure 2, (c1-c3). The recovered signal (c3) is nearly identical to f . On figure 3 one can see a real sound signal of a tiger howling, together with the signals recovered from fixed basis and adapted basis iterations. Although the final adapted basis is not the same as the one of the original signal, it still provides an improvement of 2dB with respect to a fixed spatial subdivision.

Best Bandelet Basis Compressed Sensing

Adapted Bandelet Transform

The bandelet bases dictionary was introduced by Le Pennec and Mallat [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF][START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF] to perform adaptive approximation of images with geometric singularities, such as the cartoon image in figure 4,(a). We present a simple implementation of the bandelet transform inspired from [START_REF] Peyré | Surface compression with geometric bandelets[END_REF].

A bandelet basis B λ is parameterized by λ = (Q, {θ S } S∈Q ), where Q is a quadtree segmentation of the pixels locations and θ S ∈ [0, π[ ∪ Ξ is an orientation (or the special token Ξ) defined over each square S of the segmentation, see figure 4, (a). The bandelet transform corresponding to this basis applies independently over each square S of the image either • if θ S = Ξ: a 2D isotropic wavelet transform, • if θ S = Ξ: a 1D wavelet transform along the direction defined by the angle θ S . We now detail the latter transform. The position of a pixel x = (x 1 , x 2 ) ∈ S with respect to the direction θ S is p x = sin(θ S )x 1 -cos(θ S )x 2 . The m pixels {x (i) } in S are ranked according to the 1D ordering p x (0) p x (1) . . . p x (m-1) . This ordering allows to turn the image {f ) ], see figure 4,(c). The bandelet transform of the image f inside S is defined as the 1D Haar transform of the signal f 1D , see figure 4,(d). This process is both orthogonal and easily invertible, since one only needs to compute the inverse Haar transform and pack the retrieved coefficients at the original pixels locations. Keeping only a few bandelet coefficients and setting the others to zero performs an approximation of the original image that follows the local direction θ S , see figure 4 (f).

[x]} x∈S defined over S into a 1D signal f 1D [i] = f [x (i
In order to restrict the number of tested geometries θ for a square S ∈ Q containing #S pixels, we follow [START_REF] Peyré | Surface compression with geometric bandelets[END_REF] and use the set of directions that pass through two pixels of S. The number of such directions is of the order of (#S) 2 . For this bandelet dictionary, the penalization of a basis B λ where λ = (Q, {θ S } S ) is defined as pen(λ) = #Q+ S∈Q 2 log 2 (#S), where #Q is the number of leaves in Q. A fast best basis search, described in [START_REF] Peyré | Surface compression with geometric bandelets[END_REF], allows to define a segmentation Q and a set of directions {θ S } S adapted to a given image f by minimizing (1). This process segments the image into squares S on which f is smooth, thus setting θ S = Ξ and squares containing an edge, where θ S closely matches the direction of this singularity.

Numerical Results

The geometric image depicted in figure 5, (a) is used to compare the performance of the original compressed sensing algorithm in a wavelet basis to the adaptative algorithm in a best bandelet basis. Since the wavelet basis is not The segmentation is depicted after the last iteration, together with the chosen direction θ S which closely matches the real geometry. On figure 5, (d/e/f), one can see a comparison for a natural image containing complex geometric structures such as edges, junctions and sharp line features. The best bandelet process is able to resolve these features efficiently.

Conclusion

The best basis framework presented in this paper allows to recover signals with complex structures from random measurements. This approach is successful for natural sounds and geometric images that contain a broad range of sharp transitions. Using a dictionary of bases decouples the approximation process from the redundancy needed for adaptivity and requires the design of a penalization cost on the set of bases. This lowers the computational burden and the numerical instabilities. This framework is not restricted to orthogonal bases, although it is a convenient mathematical way to ensure the compressed sensing recovery condition.

This best basis approach to sensing and recovery is also a promising avenue for interactions between biological processing, where a deterministic or chaotic process is highly probable and signal processing, where randomization has proven useful to provide universal coding strategies.

Fig. 2 .

 2 Fig. 2. (a) synthetic sound signal with 30 random cosine atoms N = 4096 ; (b) recovery using a fixed cosine basis ; (c1) first iteration of the best basis recovery algorithm, n = N/3 ; (c2) iteration s = 5 ; (c3) iteration s = 20.

Fig. 3 .

 3 Fig. 3. (a) sound signal of a tiger howling, together with the best spatial segmentation, N = 32768 ; (b) recovery using fixed local cosine basis, n = N/3 (PSNR=19.24dB) ; (c) recovery using best cosine basis, n = N/3 (PSNR=21.32dB)

Fig. 4 .

 4 Fig. 4. (a) a geometric image together with some adapted dyadic segmentation Q ; (b) a square S together with some adapted direction θS ; (c) the 1D signal f1D obtained by mapping the pixels values f (x (i) ) on a 1D axis ; (d) the 1D Haar coefficients of f1D ; (e) the 1D approximation obtained by reconstruction from the 20 largest Haar coefficients ; (f ) the corresponding square approximated in bandelet. adapted to the geometric singularities of such an image, reconstruction (b) has strong ringing artifacts. The adapted reconstruction (c) exhibits fewer such artifacts since the bandelet basis functions are elongated and follow the geometry.The segmentation is depicted after the last iteration, together with the chosen direction θ S which closely matches the real geometry. On figure5, (d/e/f), one can see a comparison for a natural image containing complex geometric structures such as edges, junctions and sharp line features. The best bandelet process is able to resolve these features efficiently.

Fig. 5 .

 5 Fig. 5. (a/d) original image ; (b/e) compressed sensing reconstruction using the wavelet basis, n = N/6 (b:PSNR=22.1dB, e:PSNR=23.2dB) ; (c/f ) reconstruction using iteration in a best bandelet basis (c: PSNR=24.3dB, f:PSNR=25.1dB).


Gabriel Peyré
email: gabriel.peyre@ceremade.dauphine.fr

Non-negative Sparse Modeling of Textures

This paper presents a statistical model for textures that uses a non-negative decomposition on a set of local atoms learned from an exemplar. This model is described by the variances and kurtosis of the marginals of the decomposition of patches in the learned dictionary. A fast sampling algorithm allows to draw a typical image from this model. The resulting texture synthesis captures the geometric features of the original exemplar. To speed up synthesis and generate structures of various sizes, a multi-scale process is used. Applications to texture synthesis, image inpainting and texture segmentation are presented.

Statistical Models for Texture Synthesis

The characterization of textures is a central topic in computer vision and graphics, mainly approached from a probabilistic point of view. Spatial domain modeling. The works of both Efros and Leung [1] and Wei and Levoy [2] pioneered a whole area of greedy approaches to texture synthesis. These methods copy pixels one by one, enforcing locally the consistence of the synthesized image with the exemplar. Recent approaches such as the method of Lefebvre and Hoppe [3] are fast, multiscale and give impressive results.

Transformed domain modeling. Julesz [START_REF] Julesz | Visual pattern discrimination[END_REF] stated simple axioms about the probabilistic characterization of textures. A texture is described as a realization of a random process characterized by the marginals of responses to a set of linear filters. Zhu, Wu and Mumford [START_REF] Zhu | Filters, random fields and maximum entropy (FRAME): Towards a unified theory for texture modeling[END_REF] setup a Gibbs energy to learn both the filters and the marginals. They use a Gibbs sampler to draw textures from this model.

A fast synthesis can be obtained by fixing the analyzing filters to be steerable wavelets as done by Heeger and Bergen [START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF]. The resulting textures are similar to those obtained by Perlin [START_REF] Perlin | An image synthesizer[END_REF]. They exhibit isotropic cloud-like structures and fail to reproduce long range anisotropic features. This is because wavelets decompositions represent sparsely point wise singularities but do not compress enough long edge features. Higher order statistics such as local correlations are used by Portilla and Simoncelli [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] to synthesize high quality textures.

Sparse image decompositions. Representing a complex image with few meaningful elements is at the core of the visual processing made by the human cortex. Atteneave [START_REF] Attneave | Some informational aspects of visual perception[END_REF] and Barlow [START_REF] Barlow | Possible principles underlying the transformation of sensory messages[END_REF] first stated that efficient high level computations should be performed over a representation of reduced complexity.

This biological processing suggests a sparse description of a patches y ∈ R N of N pixels extracted from a natural image as

y[n] = p-1 k=0 x k d k [n]
where

||x|| ℓ 0 def. = # k \ x k = 0 τ ≪ N, (1)
where x = [x 0 , . . . , x p-1 ] are the coefficients of the decomposition, and where typically p N . This simple linear model uses as prior a dictionary D = [d 0 , . . . , d p-1 ] ∈ R N ×p of atoms.

This generative process is appealing as it involves p degrees of freedom to perform statistical modeling while generating a large amount of different images. Classical works in computational harmonic analysis describe the set of candidate images y as belonging to some functional space and study the sparsity of the decomposition of y in some fixed basis or dictionary. Fourier decomposition is only suitable for smooth images and a wavelet decomposition [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] allows to model sparsely data with pointwise singularities. Images with geometrical singularities such as edges should be analyzed with more involved constructions such as the frame of curvelets [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] or a bandelet best basis [START_REF] Pennec | Bandelet Image Approximation and Compression[END_REF].

These tools however are not efficient to capture the complex geometry of textures, which might include turbulent structures or complex overlapping junctions. To sparsely represent these structures, we use a exemplar-based approach where the dictionary D is learned from a single input texture.

Learning the dictionary. Given a set of m typical patches Y = [y 0 , . . . , y m-1 ] ∈ R N ×m , one needs to compute both the dictionary D and the coefficients X = [x 0 , . . . , x m-1 ] of the decomposition Y ≈ DX that leads to the sparsity ||x j || ℓ 0 τ of each column x j of X required by equation [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF].

Olshausen and Field [START_REF] Olshausen | Emergence of simple-cell receptive-field properties by learning a sparse code for natural images[END_REF] first proposed a learning scheme to build a dictionary of atoms to represent each signal of a data set using few elements. For natural images, edge filters emerge to efficiently represent the geometry of images. Algorithms have been proposed in signal processing such as the K-SVD of Aharon et al [START_REF] Aharon | The k-svd: An algorithm for designing of overcomplete dictionaries for sparse representation[END_REF] and the tight frames construction of Tropp et al. [START_REF] Tropp | Designing structured tight frames via an alternating projection method[END_REF].

ICA and sparse dictionaries have been applied in texture modeling mainly for features extraction in classification [START_REF] Zeng | Selection of ica features for texture classification[END_REF][START_REF] Skretting | Texture classification using sparse frame based representations[END_REF]. An ICA decomposition is used as a post-processing step by Manduchi and Portilla [START_REF] Manduchi | Independent component analysis of textures[END_REF] to enhance the synthesis results of Heeger and Bergen multiscale approach [START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF].

An alternative way to perform sparse coding is to enforce positivity of both the coefficients X and the dictionary D. Non-negative factorization, as proposed by Lee and Seung [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] tends to to decompose an image into its meaningful parts, see the theoretical study of Donoho and Stodden [START_REF] Donoho | When does non-negative matrix factorization give a correct decomposition into parts[END_REF].

Contributions. In this paper, we propose a signal-processing approach to the sparse modeling of textures. It is based on the following ingredients:

Only marginal responses of a decomposition are used. It provides a simple modeling of textures ensembles in the transformed domain.

An additive generative model based on positive atoms is used. Although our method works with traditional linear decompositions, positive atoms are more localized and decompose the exemplar texture into it constitutive elements. A fast signal-processing based method is used to learn the set of atoms adapted to a given texture. This is motivated by the sparse description of textures and does not require complex solvers such as the one of [START_REF] Zhu | Filters, random fields and maximum entropy (FRAME): Towards a unified theory for texture modeling[END_REF].

A fast iterative scheme is used to sample a typical texture that matches the same marginal statistics. It does not ensure a sampling with maximum entropy distribution [START_REF] Zhu | Filters, random fields and maximum entropy (FRAME): Towards a unified theory for texture modeling[END_REF] but initializing the iterations with a random noise is good enough in typical applications, as noticed in [START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF][START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF]. The problem of learning D and X requires to minimize the reconstruction error ||Y -DX|| subject to the constraints D, X 0. This minimization is not convex on both D and X, but following [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], an alternate minimization on each matrix can be carried using iterations of the following two steps

X ab ← X ab i D ia Y ib /(DX) ib i D ia and D ia ← D ia b X ab Y ib /(DX) ib b X ab ,
which converge to a local minimum of ||Y -DX||. To enforce the dictionary elements y i to be of unit norm, the columns of D are normalized at each iteration. An explicit sparsity prior can be added into this process [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF] to enforce the constraints ||x j || ℓ 0 τ required by equation [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF]. In practice, this did not result in a noticeable enhancement for texture modeling.

Patch-based decomposition. Our sparse modeling of textures is based on the assumption of local ergodicity common to many previous works, see for instance [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF]. It assumes that the texture is a realization of a random vector whose statistics are invariant under translation.

Starting from an exemplar texture f e ∈ R N of N pixels given by the user, one works at a fixed scale w > 0 that sets the size of the typical structures present in the texture. Section 4 describes how a multiscale procedure can alleviate the issue of using a fixed scale. If n is a pixel in the image, we denote by p n = p n (f e ) the patch of w × w pixels centered at n in f e .

The ergodicity assumption leads us to model each patch p n (f e ) as being sparsely represented in some dictionary D. In practice, a set of m square patches Y = [y 0 , . . . , y m-1 ] of N = w×w pixels is extracted at random from the exemplar f e . The non negative factorization Y = DX is used to learn the dictionary from the exemplar. Figure 1 shows some examples of atoms d j of D learned from a texture.

Parametric Modeling Over a Learned Dictionary

Once a dictionary D is learned to efficiently represent the patches Y of size w × w extracted from f e , a statistical model is built using the marginal of the decomposition coefficients X such that Y = DX. Ergodicity allows to use the sets of coefficients x k = {X k,j } m-1 j=0 to estimate the marginals of the underlying probabilistic model that generates the patches of f e .

The marginal distributions are both on-sided and highly concentrated near 0. We thus keep track of the empirical variance and kurtosis of the decomposition of Y onto the dictionary D defined by

σ k (Y ) def. = σ(x k ) def. = M 2 (x k ), and κ k (Y ) def. = κ(x k ) def. = M 4 (x k ) (M 2 (x k )) 2 , (2)
where

M s (x k ) def. = 1 N m j=1 (x k [j]) s . (3)
A texture f e is characterized, at a scale w, by its adapted dictionary D, the empirical variance σ k (Y ) and kurtosis κ k (Y ). Both are computed from the decomposition of a large enough set of patches Y extracted from f e .

4 Sampling from the Positive Texture Model Texture ensembles. A dictionary D learned from the exemplar defines an equivalence relationship D ∼ between two sets of patches P = {p n } n and Q = {q n } n that shares the same statistics

P D ∼ Q ⇐⇒ ∀k, σ k (P ) = σ k (Q) and κ k (P ) = κ k (Q).
Note that these statistics are defined over a transformed domain, which means that one first has to factor the matrix P = [p 0 , p 1 , . . .] as P = DX and then extract the statistics of the rows x k , as explained in the previous section.

The set of images f of N 1 pixels whose local decompositions in D share the same marginal statistics as f e defines an ideal texture ensemble

T (f e ) def. = f \ {p n (f )} n∈P D ∼ {p n (f e )} n∈P ⊂ R N1 .
where P denotes the set of pixels. Note that although f e contains N pixels, the texture ensemble can be defined for any size N 1 . We use a model based on overlapping patches to have a translation invariant description that avoids blocking artifact.

Imposing simultaneously all the statistical constraints that define T (f e ) is a complex non-linear process due to the overlapping between patches p n (f ). We approximate this texture ensemble using the following decomposition

T (f e ) ≈ T (f e ) def. = δ∈∆ T δ (f e )
where

∆ def. = {0, . . . , w -1} 2 . ( 4
)
Each T δ (f e ) imposes constraints on a sub-set of non-overlapping patches

T δ (f e ) def. = f \ {p n (f )} n∈P δ D ∼ {p n (f e )} n∈P δ
where P δ with δ = (δ 1 , δ 2 ) is the sub-lattice of pixels n defined by

P δ = {n = (n 1 , n 2 ) \ n 1 %w = δ 1 and n 2 %w = δ 2 } ,
where % is the modulo operator.

Projection on the texture ensemble. The approximation of the texture ensemble of equation ( 4) describes T (f e ) as the intersection of convex sets T δ (f e ). Following [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], an image f can be approximately projected onto T (f e ) by iterating projections onto each T δ (f e ). Since each set T δ (f e ) involves constraints on independent patches of f , this projection can be carried by local ajustements on the decomposition of each patch p n (f ) for n ∈ P δ . We denote by π δ (f ) the approximate projection of f onto T δ (f e ), which is computed using the following steps

The set of patches P = {p n (f )} n∈P δ are gathered from f . The positive factorisation P = DX is performed using the algorithm described in section 2. We now adjust the statistics of each vector of coefficients x k = {X k,j } j representing the decomposition on a single atom d k of D. The projection on variance constraints is performed by x k ← x k σ k (f e )/σ(x k ). As done in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], enforcing the kurtosis is performed using a gradient descent of the potential x k → |κ( x k )κ k (f e )| 2 while keeping σ( x k ) constant. The updated patches P are reconstructed using the dictionary P = D X where the rows of X are the new coefficients x k for k = 0, . . . , p -1. The projection π δ (f ) is computed by rearranging the non-overlapping patches of P .

Sampling from the texture ensemble. The set T (f e ) is compact and the uniform distribution on T (f e ) thus defines the probability measure with maximum entropy. Zhu et al. [START_REF] Zhu | Filters, random fields and maximum entropy (FRAME): Towards a unified theory for texture modeling[END_REF] sample this distribution in order to synthesize textures without bias.

Rather than performing an exact sampling of the uniform distribution on T (f e ), we follow [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] and use a sampling strategy that finds a point in T (f e ) by iterating projections on each of set T δ (f e ). Starting from an initial point with high entropy such as a gaussian noise ensures that the set of generated images span T (f e ) with a minimum bias. This leads to the following synthesis algorithm:

Preprocessing: extract m random patches Y = [y 0 , . . . , y m-1 ] of width w × w from f e . Compute the positive factorization Y = DX and record the parameters σ k (f e ) and κ k (f e ) of the marginals defined in equation [START_REF] Wei | Fast texture synthesis using tree-structured vector quantization[END_REF]. Initialization: set f as a realization of a gaussian white noise on N 1 pixels. Repeat for random shift δ ∈ {0, . . . , w -1} 2 until convergence: f ← π δ (f ).

Texture Synthesis

Mono-scale synthesis. Starting from an exemplar f e of N pixels, the mono-scale texture synthesis process consists in computing an image f ∈ T (f e ) of N 1 pixel. This sampling is carried out using the iterative projection method exposed in the previous section. Note that this method generates textures of arbitrary size. Furthermore, we use cyclic boundary conditions when extracting patches {p n (f )} n∈P δ from the output texture, which results in periodic textures that tile the plane. Figure 2 shows two exemples of synthesis. The short range structures are well synthesized, but the algorithm fails to capture long range fiber-like structures.

Color texture are synthesized by applying the algorithm on each channel independently. Moving from the RGB color representation to the HSV representation improves synthesis quality since the intensity channel tends to have more distinct structures than the remaining channels. Figure 2 compares our method with the method of [START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF] and [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF]. Both the input texture f e and output f are of size 256 × 256 pixels. The multiscale histogram matching of Heeger and Bergen [START_REF] Heeger | Pyramid-Based texture analysis/synthesis[END_REF] is not able to reproduce textures with geometric features. In contrast, both the higher order model of Portaill and Simoncelli [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] and our method can synthesize textures with complex structures. Multiscale synthesis. In order to cope with the fixed scale w used in previous section, one can use a multiscale synthesis strategy. We use a fixed number of pixels w 0 but consider textures with increasing resolutions. This allows to capture first elongated low frequencies structures and then fine scale geometric details. A simple interpolation is used to switch between the various resolutions. At each scale, the synthesis algorithm manipulates only small patches of size w 0 × w 0 . This leads to the following algorithm that handles J scales. Initialization: Set j = J to be the coarser scale. Initialize the synthesis with a random noise f of N 1 /2 J × N 1 /2 J pixels.

Step 1: Set w = 2 j w 0 . Smooth the exemplar f j e = f e * h j where h j is a gaussian kernel of width 2 j pixels. Extract a set of m patches Y j from f j e . Sub-sample these squares by a factor 2 j so that vectors in Y j are of size w 0 × w 0 .

Step 2: Perform the mono-scale synthesis algorithm using patches Y j to train the dictionary and with the current f as initialization. If j = 0 then stop the algorithm. Otherwise, upsample the current synthesized texture using linear interpolation from N/2 j × N/2 j pixels to 2N/2 j × 2N/2 j pixels. Set j → j -1 and go back to step 1. In our implementation, we have used 2 scales j = 0, 1 and a base width w 0 = 8. Figure 3 compares the fixed scale synthesis and the multiscale synthesis which is able to creates elongated singularities. Figure 4 shows additional synthesis results.

Texture Inpainting

The inpainting problem consists in filling a set of missing pixels Ω in a given image f e . This problem has been approached using evolution equation derived from fluid dynamics by Bertalmio et al. [START_REF] Bertalmio | Image inpainting[END_REF] however diffusion-based approaches fail to reproduce texture patterns. Using a sparsity prior in a set of fixed bases such as Curvelets and local DCT, Fadili and Starck [START_REF] Fadili | Em algorithm for sparse representation-based image inpainting[END_REF] are able to inpaint oscillatory and elongated texture structures.

Our synthesis algorithm can be slightly modified to cope with missing data as follow.

Extract a set of m patches Y from f e that are as close as possible from Ω without intersecting it. Set as initial inpainted image f the original f e with values at random inside Ω.

Step 1: for a random shift δ, perform one step of synthesis f ← π δ (f ). The projection needs only to be performed for patches p n (f ) ∩ Ω = ∅.

Step 2: impose the known values, ∀n /

∈ Ω, f [n] ← f e [n]
. Go back to step 1. Figure 5 shows some step of this inpainting process and figure 6 shows additional results. A limitation of this method is that it works well for homogeneous textures Texture segmentation. Our model can be used to perform segmentation of a given texture f into components corresponding to patterns similar to exemplars {f 1 e , . . . , f s e }. Learned dictionaries have already been used for segmentation [START_REF] Skretting | Texture classification using sparse frame based representations[END_REF], and we recast this approach into our patch-based non-negative model. The idea is to project the texture f onto each texture ensemble T (f ℓ e ) and select locally the class ℓ that generates the least deviation from f . This leads to the following algorithm.

Learn the statistical model T (f ℓ e ) for each class ℓ. Compute the projection f ℓ of f on each ensemble T (f ℓ e ) using the algorithm of section 4. Compute the class-wise error for each pixel and smooth it with a gaussian kernel G s0 of with s 0

E ℓ [n] def. = |f [n] -f ℓ [n]| 2 and E ℓ = E ℓ * G s0 .
The smoothing removes estimation noise and reflects the prior knowledge that class boundary should be smooth curves.

Compute the segmentation into classes using

ℓ[n] = argmin ℓ E ℓ [n].
We have tested this segmentation using a set of s = 5 exemplar textures. The input image f of 256 × 256 pixels is a patchwork of five textures extracted from the upper left corner of the original images f ℓ e of 512 × 512 pixels. The exemplars f ℓ e are extracted from the lower right corner of f ℓ e . Figure 7 shows the segmentation process.

Conclusion

We have proposed a statistical model for textures built out of marginal distributions of a positive decomposition. The statistical model is parametric since we use only low order moments of the distribution. This is permitted thanks to the sparsity provided by a learned dictionary. Such a positive dictionary captures with few atoms the structures of the textures. This simple model allows to perform multiscale texture synthesis and can be fitted into various applications such as texture inpainting or texture segmentation.

An important parameter of our model is the redundancy factor p/N . Redundancy brings invariances to various factors such as translations or local illumination changes. These parameters are however hard to control. Representations that can explicitly capture this invariances include bilinear decompositions [START_REF] Grimes | Bilinear sparse coding for invariant vision[END_REF] that could improve our model.

2

  Non-negative Atoms for Texture DecompositionNon-negative matrix factorization. The process of learning a dictionary D of atoms D = [d 0 , . . . , d p-1 ] to represent accurately a set Y = [y 0 , . . . , y m-1 ] of exemplars is equivalent to performing a factorization Y = DX. This problem is underconstrained and the non-negative factorization[START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] enforces positivity of both X and D to constrained the learning problem.

Fig. 1 .

 1 Fig. 1. Example of dictionary for patches of 12 × 12 pixels. We used m = 20w 2 random patches for the learning stage.

Fig. 2 .

 2 Fig. 2. (a) Original texture. (b) Textures synthesized with the method of Heeger and Bergen [6] (c) Textures synthesized with the method of Portilla and Simoncelli [8]. (d) Textures synthesized with our method.

Fig. 3 .

 3 Fig. 3. (a) Original texture fe. (b) Texture synthesized with the mono-scale procedure with patches of width w = 8. (c) Texture synthesized at scale 2 j = 2 with w = 16. (d) Texture synthesized at scale 2 j = 1 with w = 8.

Fig. 4 .

 4 Fig. 4. Examples of multiscale texture synthesis.

  and the inpainting tends to give poor results if Ω intersects a broad range of structures.

Fig. 5 .

 5 Fig. 5. (a) Texture to inpaint, the missing region Ω is depicted in black. (b,c,d) Evolution of the inpainting for step 1,2,4. (e) Final result.

Fig. 6 .

 6 Fig. 6. Examples of inpainting.

Fig. 7 .

 7 Fig. 7. (a) Original texture f . (b) Projected texture f1 of f onto T (f 1 e ). Note how the upper left corner is well preserved. (c) Projected texture f2. (d) Projected texture f3. (e) Ground trust segmentation. (f ) Segmentation ℓ[n] computed with s0 = 3 pixels. (g) Segmentation computed with s0 = 6 pixels.


Mehmet Süzen
email: suzen@acm.org

Cooperative Phenomenon: Statistical Physics Origins of Connectionist Learning

Keywords: Neural Networks, Statistical Physics, Neurophysiology, Spin-glasses, Lattice Models, Dynamics, Learning, Statistics of Cooperative Phenomenon, Lenz-Ising Architectures

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Understanding natural or artificial phenomenon in the language of discrete mathematics is probably one of the most powerful toolbox scientist use [START_REF] Kenneth H Rosen | Handbook of Discrete and Combinatorial Mathematics[END_REF]. Large portion of computer science and statistical physics deals with such finite structures. One of the most prominent successful usage of such approach was Lenz and Ising's work on modelling ferromagnetic materials [START_REF] Lenz | Beitrag zum Verstländnis der Magnetischen Erscheinungen in Festen Körpern[END_REF][START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF][START_REF] Ising | The fate of Ernst Ising and the fate of his model[END_REF][START_REF] David | A guide to Monte Carlo Simulations in Statistical Physics[END_REF] and neural networks as a model to biological neuronal structures [START_REF] Mcculloch | A Logical Calculus of the Ideas Imminent in Nervous Activity[END_REF][START_REF] Cole | Representation of Events in Nerve Nets and Finite Automata[END_REF][START_REF] Little | The Existence of Persistent States in the Brain[END_REF].

The analogy between two areas of distinct research have been pointed out by many researchers [START_REF] Peretto | Collective Properties of Neural Networks: a Statistical Physics Approach[END_REF][START_REF] Van Hemmen | Spin-glass Models of a Neural Network[END_REF][START_REF] Sompolinsky | Statistical Mechanics of Neural Networks[END_REF][START_REF] Sherrington | Neural Networks: the Spin Glass Approach[END_REF][START_REF] Bahri | Statistical Mechanics of Deep Learning[END_REF]. However, the discourse and evolution of these approaches were kept as two distinct research fields and many innovative approaches rediscovered under different names.

Cooperative phenomenon

Statistical definition of cooperative phenomenon pioneered by Wannier and Kremer [START_REF] Gregory | The Statistical Problem in Cooperative Phenomena[END_REF][START_REF] Hendrik | Statistics of the two-dimensional ferromagnet[END_REF][START_REF] Hendrik | Statistics of the two-dimensional ferromagnet[END_REF]. Even though their technical work focused on extension of Ising model to 2D with cyclic boundary condition and introduction of exact solutions with matrix algebra, they were the first to document the potential of how Lenz-Ising model actually represent a more generic system than marely model to ferromagnets, namely anything falls under cooperative phenomenon can be addressed with Lenz-Ising type model, summarized in Definition 1.

Definition 1. Cooperative phenomenon of Wannier type [START_REF] Gregory | The Statistical Problem in Cooperative Phenomena[END_REF] Set of N discrete units, U , identified with a function s i , i=1,..,N forms a collection or assembly. The function that identifies the units is a mapping s i : R → R. A statistic S applied on U is called cooperative phenomenon of Wannier type W .

A statistic S can be any mapping or set of operations on the assembly of units U . For example inducing ordering on the assembly of units and summation over s i values, would corresponds to non-interacting magnetic system with unit external field or non-connected set of neurons capacity of inhibition or exhibition. However, amazingly, Definition 1 is so generic that Rosenblatt's perceptron [START_REF] Van Der Malsburg | Frank Rosenblatt: principles of neurodynamics: perceptrons and the theory of brain mechanisms[END_REF], current deep learning systems [START_REF] Schmidhuber | Deep learning in Neural Networks: An overview[END_REF] and complex networks [START_REF] Duncan | Collective dynamics of 'small-world'networks[END_REF] falls into this category as well.

The originality of Cooperative phenomenon of Wannier type comes on a secondary concept, so called event propagation as given in Definition 2.

Definition 2. Event propagation [START_REF] Gregory | The Statistical Problem in Cooperative Phenomena[END_REF] An event is defined as a snapshot of cooperative phenomenon of Wannier type W . If an event takes place of one unit of assembly U , the same event will be favored by other units, this is expressed as event propagation between two disjoint set of units E (u 1 , u 2 ), and u 1 ∩ u 2 = ∅ and u 1 , u 2 ∈ U and with an additional statistic S is defined.

The parallels between Wannier's event propagations are remarkably the same as of neural network formalism defined by McCulloch-Pitts-Kleene [START_REF] Mcculloch | A Logical Calculus of the Ideas Imminent in Nervous Activity[END_REF][START_REF] Cole | Representation of Events in Nerve Nets and Finite Automata[END_REF], not only conceptually but matematical treatment is identical and originates from Lenz-Ising model's treatment of discrete units. As we mentioned, this goes beyond doubt not a simple analogy but forms a generic framework as envisioned by Wannier. The similarity between ferromagnetic systems and neural networks is probably first documented directly by Little [START_REF] Little | The Existence of Persistent States in the Brain[END_REF]: Spin states of magnetic spins corresponds to firing state of a neuron. Unfortunately, Little only see it as simple analogy, and missed the opportunity provided by Wannier as a generic natural phenomenon of cooperation.

The conceptual similarity and inference on Wannier's event propagation appears to be quite close to Hebb's learning [START_REF] Olding | The Organization of Behavior: a Neuropsychological Theory[END_REF] and gives natural justification for backpropagation for multilayered networks. History of backpropagation is exhaustively studied elsewhere [START_REF] Schmidhuber | Deep learning in Neural Networks: An overview[END_REF].

Lenz-Ising architectures: Ferromagnets to Nerve Nets

As we establised two basic definitions of cooperative phenomenon, we can now define a generic setting of Lenz-Ising model that captures both physics literature that extensively used this in so called spin-glasses research and for neural networks. A guiding principle will be based on Wannier's definition of cooperative phenomenon.

Definition 3. Lenz-Ising Architectures

Given Wannier type cooperative phenomenon W , imposing constrains on the discrete units, U c that they should be spatially correlated on the edges E of an arbirary graph G (E, V ) with ordering and with vertices V of the arbitrary graph carring coupling weight between connected two units with biases. Set of event propagations E c defined on the cooperative phenomeon can induce dynamics on defining vertice weights, or vice versa. ILAs are defined as statistic S applied to U c with propagations E c . Lenz-Ising Architectures (ILAs) should not be confused with graph neural networks as it does not model data structures. It could be seen as subset of graph dynamical systems in some sense but formal connections should be established elsewhere. However, primary characteristic of ILAs are that it is conceptual and mathematical representation of spin-glass systems (including Lenz-Ising, Anderson, Sherrington-Kirkpatrick, Potts systems) and neural networks (including recurrent and convolutional networks) under the same umbrella.

Learning representations inherent in Metropolis-Glauber dynamics

The primary originality in any neural network research papers lies in so called learning representation from data and generalisation. However, it isn't obvious to the that community that actually spin-glasses are capable of learning representations inherently by induced dynamics such as Metropolis or Glauber dynamics by construction, as an inverse problem.

In physics literature this appears as finding a solution to the problem of how to express free energy and minimisation of this with respect to weights or coupling coefficients, This is noting but a learning represenations. Usually a simulation approach is taken as a route, for example Monte Carlo techniques [START_REF] David | A guide to Monte Carlo Simulations in Statistical Physics[END_REF][START_REF] Süzen | Effective ergodicity in single-spin-flip dynamics[END_REF][START_REF] Süzen | Anomalous diffusion in convergence to effective ergodicity[END_REF] via Metropolis or Glauber dynamics. The intimate connection between concepts of ergodicity and learning in deep learning is recently shown [START_REF] Bahri | Statistical Mechanics of Deep Learning[END_REF][START_REF] Süzen | Spectral ergodicity in deep learning architectures via surrogate random matrices[END_REF][START_REF] Süzen | Periodic Spectral Ergodicity: A Complexity Measure for Deep Neural Networks and Neural Architecture Search[END_REF] in this context.

As we argued earlier the generic definition provided by Wannier on cooperative phenomenon and ILAs; there is an intimate connection with learning and so called solving spin-glasses that usually boils down to computing free energies as mentioned. And a link between two distinct fields, computing backpropagation and free energies are natural candidates to establish equivalance relations.

Conclusions and Outlook

Apart from honouring physicists Lenz and Ising, based on understanding of cooperative phenomenon's origins, naming the research outpus from of spin-glasses and neural networks under an umbrella term Lenz-Ising architectures (ILAs) is historically accurate and technically a resonable naming scheme under the overwhelming evidence given in the literature. This is akin to naming current computers with von Neumann architectures. This forms the origins of connectionist learning from statistical physics, where this approach currently enjoying vast engineering success today.

The rich connection between two areas in computer science and statistical physics should be celebrated. For more fruitful collaborations, both literatures, embracing large statistics literature as well, should converge much more closely. This would help communities to avoid awkward situations of reinventing the wheel again and hindering recognition of the work done by physicists decades earlies, i.e., Ising and Lenz. this manuscript. Author has no competing or other kind of conflict of interests. This work is produced solely with the aim of scholarly work and does not have any personal nature at all. Manuscript is dedicated in memory of Ernst Ising for his contribution to physics of ferromagnetic materials, now seems to have far more implications.

Acknowledgements and Notes

Author is grateful to Jürgen Schmidhuber for correspondence, pointing out the work of S.C.Kleene on the dynamics of neural networks and encouraging the inceptions of


Oanh Chau
email: oanh.chau@univ-reunion.fr

Adrien Petrov
email: apetrov@math.univ-lyon1.fr

Arnaud Heibig
email: arnaud.heibig@insa-lyon1

A class of dynamic unilateral contact problems with sub-differential friction law

Keywords: Mathematics Subject Classification (2000). 74M15, 74M10, 74F05, 74H20, 74H25, 34G25 time depending thermo-visco-elasticity, unilateral contact, sub-differential friction law, non clamped condition, evolution variational inequality

We study a class of dynamic unilateral contact problems with subdifferential friction law, and thermal effects, for time depending long memory viscoelastic materials, with or without the clamped condition. We describe the mechanical problem, derive its variational formulation, and after specifying the assumptions on the data and operators, we prove an existence and uniqueness of weak solution on displacement and temperature fields.

Introduction

Since the dawn of time, situations of contact between deformable bodies abound in everyday life and engineering applications, where the numerous forces acting on the system may lead to the appearance of microcracks, and deteriorate the mechanical equipments. It is then important to understand the complexity of the contact phenomena, in order to guarantee the safety of the mechanism of the system. There-fore considerable efforts have been achieved in modeling, mathematical analysis and numerical simulations, within the weak distributional formulation framework, expressed in terms of evolutional variational inequalities and hemivariational inequalities. The literature dedicated to this field is increasing day by day. The state of the art can be found in the masterpieces [START_REF] Duvaut | Les Inéquations en Mécanique et en Physique[END_REF], [START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Applications[END_REF], [START_REF] Panagiotopoulos | Hemivariational Inequalities, applications in Mechanics and Engineering, Contributions to Nonlinear Functional Analysis[END_REF].

This work is a continuation of the paper in [START_REF] Chau | Solvability for a thermodynamic unilateral contact problem with friction and damage[END_REF], where the authors studied a frictional unilateral contact problem undergoing the thermal expansion and damage, for viscoelastic solid. The friction is modeled by the Coulomb's dry friction law.

Here we investigate a thermal contact with friction, for time depending long memory visco-elastic materials, with or without the clamped condition. The contact conditions are unilateral. The friction is modeled by an inequality of sub-differential type.

The paper is organized as follows. In Section 2 we describe the mechanical problem, its corresponding variational formulation, and then we claim the main existence and uniqueness result under specific assumptions, that we prove in Section 3.

The contact problem

In this section we describe the mechanical problem, list the assumptions on the data and derive the corresponding variational formulation. Then we state an existence and uniqueness result on displacement field and temperature, which we will prove in the next section. The physical setting is as follows. A visco-elastic body occupies a domain Ω in R d (d = 1, d = 2 or d = 3) with a Lipschitz boundary Γ that is partionned into three disjoint measurable parts, Γ 1 , Γ 2 and Γ 3 . We denote by ν the unit outward normal on Γ . Let [0, T ] be the time interval of interest, where T > 0. The body is clamped on Γ 1 × (0, T ) and therefore the displacement field vanishes there. Here we suppose that meas(Γ 1 ) = 0 or meas(Γ 1 ) > 0, which means that Γ 1 may be an empty set or reduced to a finite set of points. We assume that a volume force of density f 0 acts in Ω × (0, T ) and that surface tractions of density f 2 act on Γ 2 × (0, T ). The body may come in contact with an obstacle, the foundation, over the contact surface Γ 3 . The contact is unilateral involving here both the normal displacement, the normal velocity and the normal component of the Cauchy stresss vector. The model of the friction is specified by a general sub-differential condition, where thermal effects may occur in the frictional contact with the basis. We are interested in the dynamic evolution of the body. Let us recall now some classical notations, see e.g. [START_REF] Duvaut | Les Inéquations en Mécanique et en Physique[END_REF] for further details. We denote by S d the space of second order symmetric tensors on R d , while " • " and | • | will represent the inner product and the Euclidean norm on S d and R d . Everywhere in the sequel the indices i and j run from 1 to d, summation over repeated indices is implied and the index that follows a comma represents the partial derivative with respect to the corresponding component of the independent variable. We also use the following notation:

H = L 2 (Ω ) d , H = { σ = (σ i j ) | σ i j = σ ji ∈ L 2 (Ω ), 1 ≤ i, j ≤ d}, H 1 = { u ∈ H | ε(u) ∈ H }, H 1 = { σ ∈ H | Div σ ∈ H }.
Here ε : H 1 -→ H and Div : H 1 -→ H are the deformation and the divergence operators, respectively, defined by :

ε(u) = (ε i j (u)), ε i j (u) = 1 2 (u i, j + u j,i ), Div σ = (σ i j, j ).
The spaces H, H , H 1 and H 1 are real Hilbert spaces endowed with the canonical inner products given by :

(u, v) H = Ω u i v i dx, (σ , τ) H = Ω σ i j τ i j dx, (u, v) H 1 = (u, v) H + (ε(u), ε(v)) H , (σ , τ) H 1 = (σ , τ) H + (Div σ , Div τ) H .
Recall that D(Ω ) denotes the set of infinitely differentiable real functions with compact support in Ω ; and W m,p (Ω ), H m (Ω ) := W m,2 (Ω ), m ∈ N, 1 ≤ p ≤ +∞ for the classical real Sobolev spaces; L p (U; X) the classical L p spaces defined on U with values in X.

To continue, the mechanical problem is then formulated as follows.

Problem Q : Find a displacement field u : (0, T ) × Ω -→ R d and a stress field σ : (0, T ) × Ω -→ S d and a temperature field θ : (0, T ) × Ω -→ R + such that for a.e. t ∈ (0, T ):

σ (t) = A (t)ε(u(t)) + G (t)ε(u(t)) + t 0 B(t -s) ε(u(s)) ds + C e (t, θ (t)) in Ω (1.1) ü(t) = Div σ (t) + f 0 (t) in Ω (1.2) u(t) = 0 on Γ 1 (1.3) σ (t)ν = f 2 (t) on Γ 2 (1.4) u ν (t) ≤ 0, σ ν (t) ≤ 0, σ ν (t) • uν (t) = 0 on Γ 3 , (1.5)
ϕ c (t, w τ ) -ϕ c (t, uτ (t)) ≥ -σ τ (t) • (w τ -uτ (t)) ∀w ∈ V on Γ 3 (1.6) θ (t) -div(K c (t, ∇θ (t))) = D e (t, u(t)) + q(t) in Ω , (1.7)
-K c (t, x, ∇θ (t, x)) ν := ϕ thermal (t, x, θ (t, x)) a.e. x ∈ Γ 3 , (1.8)
θ (t) = 0 on Γ 1 ∪ Γ 2 (1.9) θ (0) = θ 0 in Ω (1.10) u(0) = u 0 , u(0) = v 0 in Ω (1.11)
Here, (1.1) is the Kelving Voigt's time-dependent long memory thermo-visco-elastic constitutive law of the body, where σ represents the stress tensor; A denotes the viscosity operator depending on the velocity of infinitesimal deformations ε(u), with the notation : for τ ∈ S d , A (t)τ = A (t, •, τ) some function defined on Ω ; here a dot above a quantity represents the derivative of the quantity with respect to the time variable; G is the elastic operator depending on the linearized strain tensor ε(u) of infinitesimal deformations, with G (t)τ = G (t, •, τ) which is defined on Ω . The term B(t)τ = B(t, •, τ) represents the so called relaxation tensor which is time-depending on the linearized strain tensor and is defined on Ω . Recall that the visco-elastic short memory corresponds to the case B ≡ 0. The last tensor C e (t, θ ) := C e (t, •, θ ) denotes the thermal expansion tensor depending on time and on the temperature, defined on Ω . For example,

C e (t, θ ) := -θ C exp (t) in Ω , where C exp (t) := (c i j (t, •))
is some time-depending expansion tensor, defined on Ω . In (1.2) is the dynamic equation of motion where the mass density ρ ≡ 1. The equation in (1.3) is the clamped condition. In (1.4) is the traction condition, where σ ν represents the Cauchy stress vector.

On the contact surface, the general relations in (1.5) are the so-called Signorini's boundary conditions (see [START_REF] Signorini | Questioni di elasticità non linearizzata e semilinearizzata[END_REF]), involving here both the normal displacement u ν , the normal velocity uν and the normal component of the Cauchy stresss vector σ ν . Recall that the condition u ν ≤ 0 represents a non penetration of the surface asperities into the obstacle : the inequality u ν < 0 means that there is no contact of the surface with the obstacle, and the equality u ν = 0 denotes contact. The last condition stipulates that the product σ ν uν vanishes almost everywhere on (0, T ) × Γ 3 . It's physical meaning can be seen as a natural consequence of the non penetrability, under some regularity assumption. To show that, let fix

(t, x) ∈ (0, T ) × Γ 3 . Case 1 : u ν (t, x) < 0. There is no contact, then σ ν (t, x) = 0 and σ ν uν = 0 at (t, x). Case 2 : u ν (t, x) = 0. Consider three sub-cases : (i) uν (t, x) > 0; (ii) uν (t, x) = 0; (iii) uν (t, x) < 0.
In the case (i), as u ν (t+h,x)-u ν (t,x) h → uν (t, x), h → 0, h > 0; then u ν (t + h, x) > 0 for h > 0 small enought, which contradicts the non penetration. In the case (ii), we have σ ν uν = 0 at (t, x). In the case (iii), as in (i) we deduce that u ν (t + h, x) < 0 for any h > 0 small enought, thus there is no contact and σ ν (t + h, x) = 0 for any h > 0 small enought, and then the continuity regularity implies that σ ν (t, x) = 0 and σ ν uν = 0 at (t, x).

The friction on the contact surface is modeled by the equation in (1.6), which is a sub-differential boundary condition, in the classical framework of convex analysis. Here, V denotes the space of admissible displacements which will be detailed in the following; σ τ represents the tangential component of the Cauchy stress vector; and ϕ c : (0, T ) × Γ 3 × R d -→ R is a given sub-differential friction contact function. Various situations may be modelled by such a condition, see below at the end of this Section. The differential equation (1.7) describes the evolution of the temperature field, where K c (t, ∇θ ) := K c (t, •, ∇θ ) is some nonlinear thermal conductivity function defined on Ω , depending on time and on the temperature gradient ∇θ . For example, denote by

K c (t, •) := (k i j (t, •))
the thermal conductivity tensor defined on Ω , we could consider

K c (t, •, ∇θ ) = K c (t, •) ∇θ .
In the second member, D e (t, u(t)) := D e (t, •, u(t)) represents some deformationviscosity heat nonlinear function defied on Ω and depending on the displacement velocity, whereas q(t) denotes the density of volume heat sources. For example,

D e (t, u(t)) = -C exp (t) : u(t) = -c i j (t, •) ∂ ui ∂ x j (t).
The associated temperature boundary condition is given by (1.8) and (1.9), where ϕ thermal are some thermal boundary function defined on (0, T ) × Γ 3 × R. Here

ϕ thermal (t, x, r) := [ϕ thermal (t, x, •)] (r), ∀(t, x, r) ∈ (0, T ) × Γ 3 × R
denotes the derivative on the third variable of ϕ thermal .

Taking the previous example for K c , we have

K c (t, x, ∇θ ) ν = k i j (t, x) ∂ θ ∂ x j ν i .
Let consider the following standard example

ϕ thermal (t, x, r) := 1 2 k e (t, x)(r -θ R (t, x)) 2 , ∀(t, x, r) ∈ (0, T ) × Γ 3 × R, (1.12)
where θ R is the temperature of the foundation, and k e is the heat exchange coefficient between the body and the obstacle. We obtain

ϕ thermal (t, x, r) = k e (t, x) (r -θ R (t, x)), (t, x, r) ∈ (0, T ) × Γ 3 × R.
Finally in (1.10) and (1.11), θ 0 , u 0 , v 0 represent the initial temperature, displacement and velocity respectively.

One may remark that since ϕ c is assumed real-valued, then unilateral contact, defined by indicator functions taking infinite values, is excluded. So the body is in fixed contact with the foundation of the body according to a friction law. This is consistent with the linear heat conduction modeled in (1.7). We insist that the new feature here is that we may have the absence of the usual claimed condition in the case where meas(Γ 1 ) = 0. However, there is coerciveness with regard to the temperature by (1.8). To derive the variational formulation of the mechanical problems (1.1)-(1.11) we need additional notations. Thus, let consider the space V of admissible displacements defined by the following closed subspace of

H 1 V = { w ∈ H 1 | w = 0 on Γ 1 and w ν ≤ 0 }.
On V we consider the inner product given by

(u, v) V = (ε(u), ε(v)) H + (u, v) H ∀ u, v ∈ V,
and let • V be the associated norm, i.e.

v 2 V = ε(v) 2 H + v 2 H ∀ v ∈ V.
It follows that • H 1 and • V are equivalent norms on V and therefore (V, • V ) is a real Hilbert space. Moreover, by the Sobolev's trace theorem, we have a constant C 0 > 0 depending only on Ω , and Γ 3 such that

v L 2 (Γ 3 ) ≤ C 0 v V ∀ v ∈ V.
Consider then the following spaces for the temperature field:

E = {η ∈ H 1 (Ω ), η = 0 on Γ 1 ∪ Γ 2 }; F = L 2 (Ω ).
The spaces E and F, endowed with their respective canonical inner product, are Hilbert spaces.

Identifying then H and F with their own duals, we obtain two Gelfand evolution triples (see e.g. [START_REF] Zeidler | Nonlinear Functional Analysis and its Applications[END_REF] II/A p. 416):

V ⊂ H ≡ H ⊂ V , E ⊂ F ≡ F ⊂ E
where the inclusions are continuous and dense. Finally, we use the notation •, • V ×V and •, • E ×E to represent the duality pairing between V and V , and respectively between E and E, which means :

u, v V ×V = u, v H , ∀ u ∈ H, ∀ v ∈ V. and η, ξ E ×E = η, ξ F , ∀ η ∈ F, ∀ ξ ∈ E.
In the study of the mechanical problem (1.1)-(1.11), we assume that the viscosity operator A :

(0, T ) × Ω × S d -→ S d , (t, x, τ) -→ A (t, x, τ) satisfies                              (i) A (•, •, τ) is measurable on (0, T ) × Ω , ∀τ ∈ S d ;
(ii) A (t, x, •) is continuous on S d for a.e. (t, x) ∈ (0, T ) × Ω ;

(iii) there exists m A > 0 such that

(A (t, x, τ 1 ) -A (t, x, τ 2 )) • (τ 1 -τ 2 ) ≥ m A |τ 1 -τ 2 | 2 , ∀τ 1 , τ 2 ∈ S d , for a.e. (t, x) ∈ (0, T ) × Ω ; (iv) there exists c A 0 ∈ L 2 ((0, T ) × Ω ; R + ), c A 1 > 0 such that |A (t, x, τ)| ≤ c A 0 (t, x) + c A 1 |τ|, ∀τ ∈ S d , for a.e. (t, x) ∈ (0, T ) × Ω . (1.13)
In this paper for every t ∈ (0, T ), τ ∈ S d we denote by A (t) = A (t, •, •) a functional which is defined on Ω × S d and A (t) τ = A (t, •, τ) some function defined on Ω . The elasticity operator G : (0, T ) × Ω × S d -→ S d satisfies :

               (i) there exists L G > 0 such that |G (t, x, ε 1 ) -G (t, x, ε 2 )| ≤ L G |ε 1 -ε 2 | ∀ε 1 , ε 2 ∈ S d , a.e. (t, x) ∈ (0, T ) × Ω ; (ii) G (•, •, ε) is Lebesgue measurable on (0, T ) × Ω , ∀ε ∈ S d ; (iii) the mapping G (•, •, 0) ∈ H . (1.14)
We put again G (t)τ = G (t, •, τ) some function defined on Ω for every t ∈ (0, T ), τ ∈ S d . The relaxation tensor B : where we denote by B(t)τ = B(t, •, τ) which is defined on Ω for every t ∈ (0, T ), τ ∈ S d . We suppose the body forces and surface tractions satisfy

(0, T ) × Ω × S d -→ S d , (t, x, τ) -→ (B i jkh (t, x) τ kh ) sat- isfies      (i) B i jkh ∈ L ∞ ((0, T ) × Ω ); (ii) B(t)σ • τ = σ • B(t)τ ∀σ , τ ∈ S d
f 0 ∈ L 2 (0, T ; H), f 2 ∈ L 2 (0, T ; L 2 (Γ 2 ) d ) (1.16)
On the contact surface, the following frictional contact function

ψ c : (0, T ) ×V -→ R, ψ c (t, w) := Γ 3 ϕ c (t, w) da, ∀(t, w) ∈ (0, T ) ×V, verifies        (i) t ∈ (0, T ) -→ ψ c (t, w) is Lebesgue measurable ∀w ∈ V ; (ii) |ψ c (t, w)| ≤ c(t) + d w V , ∀w ∈ V, a.e. t ∈ (0, T ); (iii) ψ c (t, •) is convex on V a.e. t ∈ (0, T ), (1.17)
where d > 0 is some constante and c ∈ L 2 (0, T ; R + ).

The thermal expansion tensor C e : (0,

T ) × Ω × R -→ S d verifies                      (i) C e (•, •, ϑ ) is measurable on (0, T ) × Ω , ∀ϑ ∈ R;
(ii) there exists L e > 0 such that |C e (t, x, ϑ 1 ) -C e (t, x,

ϑ 2 )| ≤ L e |ϑ 1 -ϑ 2 | ∀ϑ 1 , ϑ 2 ∈ R, a.e. (t, x) ∈ (0, T ) × Ω ; (iii) there exists c C e 0 ∈ L ∞ ((0, T ) × Ω ; R + ), c C e 1 ≥ 0 such that |C e (t, x, ϑ )| ≤ c C e 0 (t, x) + c C e 1 |ϑ |, ∀ϑ ∈ R, for a.e. (t, x) ∈ (0, T ) × Ω .
(1.18) Here we recall the notation C e (t, ϑ ) = C e (t, •, ϑ ) some function defined on Ω , for all t ∈ (0, T ) and ϑ ∈ R.

The nonlinear function K

c : (0, T ) × Ω × R d -→ R satisfies :                                  (i) K c (•, •, ξ ) is measurable on (0, T ) × Ω , ∀ξ ∈ R d ; (ii) K c (t, x, •) is continuous on R d , a.e. (t, x) ∈ (0, T ) × Ω ; (iii) there exists c K c 0 ∈ L 2 ((0, T ) × Ω ; R + ), c K c 1 ≥ 0, such that |K c (t, x, ξ )| ≤ c K c 0 (t, x) + c K c 1 |ξ |, ∀ξ ∈ R d , a.e. (t, x) ∈ (0, T ) × Ω ; (iv) there exists m K c > 0 such that (K c (t, x, ξ 1 ) -K c (t, x, ξ 2 )) • (ξ 1 -ξ 2 ) ≥ m K c |ξ 1 -ξ 2 | 2 , ∀ξ 1 , ξ 2 ∈ R d , a.e. (t, x) ∈ (0, T ) × Ω ; (v) there exists n K c > 0 such that K c (t, x, ξ ) • ξ ≥ n K c |ξ | 2 , ∀ξ ∈ R d , a.e. (t, x) ∈ (0, T ) × Ω . (1.19) We suppose that the deformation-viscosity heat D e : (0, T ) × Ω × R d -→ R satisfies :            (i) D e (•, •, v) is measurable on (0, T ) × Ω , ∀v ∈ R d ; (ii) there exists L D e > 0 such that |D e (t, x, v 1 ) -D e (t, x, v 2 )| ≤ L D e |v 1 -v 2 |, ∀v 1 , v 2 ∈ R d , a.e. (t, x) ∈ (0, T ) × Ω .
(1.20)

We assume for the heat sources density, that

q ∈ L 2 (0, T ; L 2 (Ω )) (1.21)
The nonlinear function

ϕ thermal : (0, T ) × Γ c × R -→ R verifies :                  (i) ϕ thermal (•, •, r) is measurable on (0, T ) × Γ c , ∀r ∈ R;
(ii) ϕ thermal (t, x, •) is convex derivable on R for a.e. (t, x) ∈ (0, T ) × Γ c ;

(iii) there exists c

ϕ 0 ∈ L 2 ((0, T ) × Γ c ; R + ), c ϕ 1 ≥ 0, such that |ϕ thermal (t, x, r)| ≤ c ϕ 0 (t, x) + c ϕ 1 |r|, ∀r ∈ R, a.e. (t, x) ∈ (0, T ) × Γ c . (1.22)
We notice that these assumptions are verified for the example (1.12). Finally we assume that the initial data satisfy the conditions

u 0 ∈ H, v 0 ∈ V, θ 0 ∈ E. (1.23)
To continue, using Green's formula, we obtain the variational formulation of the mechanical problem Q in abstract form as follows.

Problem QV : Find u : [0, T ] → V , θ : [0, T ] → E satisfying a.e. t ∈ (0, T ):        ü(t) + A(t) u(t) + B(t) u(t) +C(t) θ (t), w -u(t) V ×V +( t 0 B(t -s) ε(u(s)) ds, ε(w) -ε(u(t))) H + ψ c (t, w τ ) -ψ c (t, uτ (t)) ≥ f (t), w -u(t) V ×V ∀w ∈ V. (1.24) θ (t), η E ×E + K(t) θ (t), η E ×E + ψ thermal (t) θ (t), η E ×E = R(t)u(t), η E ×E + Q(t), η E ×E , ∀ η ∈ E. (1.25) u(0) = u 0 , u(0) = v 0 , θ (0) = θ 0 . (1.26)
Here, the operators and functions A(t), B(t

) : V -→ V , C(t) : E -→ V , f : [0, T ] -→ V , K(t) : E -→ E , ψ thermal (t) : E -→ E , R(t) : V -→ E and Q : [0, T ] -→ E are defined by ∀v ∈ V , ∀w ∈ V , ∀ζ ∈ E, ∀η ∈ E, a.e. t ∈ (0, T ): A(t) v, w V ×V = (A (t) ε(v), ε(w)) H ; (1.27) B(t) v, w V ×V = (G (t) ε(v), ε(w)) H ; (1.28) C(t)ζ , w V ×V = (C e (t, ζ ), ε(w)) H ; (1.29) f (t), w V ×V = ( f 0 (t), w) H + ( f 2 (t), w) (L 2 (Γ 2 )) d ; (1.30)
K(t) ζ , η E ×E = Ω K c (t, ∇ζ ) • ∇η dx; (1.31)
ψ thermal (t) ζ , η E ×E = Γ 3 ϕ thermal (t, x, ζ (x)) η(x) da(x); (1.32)
R(t) v, η E ×E = Ω D e (t, v) η dx; (1.33) Q(t), η E ×E = Ω q(t) η dx. (1.34)
Our main existence and uniqueness result is stated as follows, that we prove in the next Section.

Theorem 1. Assume that (1.13)-(1.23) hold, then there exists an unique solution {u, θ } to the problem QV with the regularity :

u ∈ W 1,2 (0, T ;V ) ∩W 2,2 (0, T ;V ) ∩C 1 (0, T ; H) θ ∈ L 2 (0, T ; E) ∩W 1,2 (0, T ; E ) ∩C(0, T ; F). (1.35)
Before proving the main theorem, we present here some examples with subdifferential friction laws of the form (1.6), see e.g. the monograph [START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Applications[END_REF] or the Habilitation thesis [START_REF] Chau | Quelques problèmes d'évolution en mécanique de contact et en biochimie[END_REF] p. 117.

Example 1. Contact with Tresca's friction law. This contact condition can be found in [START_REF] Duvaut | Les Inéquations en Mécanique et en Physique[END_REF][START_REF] Panagiotopoulos | Hemivariational Inequalities, applications in Mechanics and Engineering, Contributions to Nonlinear Functional Analysis[END_REF]. It is in the form of the following boundary condition:

   u ν = 0, |σ τ | ≤ g, |σ τ | < g =⇒ uτ = 0, on Γ 3 × (0, T ). |σ τ | = g =⇒ uτ = -λ σ τ , λ ≥ 0 (1.36)
Here g ∈ L ∞ (Γ 3 ; R + ) represents the friction bound, i.e., the magnitude of the limiting friction traction at which slip begins. The contact is assumed to be bilateral, i.e., there is no loss of contact during the process. We choose here ϕ c (v) = g |v|, ∀v ∈ V.

Example 2. Contact with viscoelastic friction condition. We consider the problems with the boundary conditions

u ν = 0, σ τ = -µ|u τ | p-1 uτ on Γ 3 × (0, T ), (1.37)
where µ ∈ L ∞ (Γ 3 ; R + ) is the coefficient of friction and 0 < p ≤ 1. Here, the tangential shear is proportional to the power p of the tangential speed, which is the case when the contact surface is lubricated with a thin layer of non-Newtonian fluid.

Here we choose

ϕ c (v) = µ p + 1 |v| p+1 , ∀v ∈ V.

Proof of Theorem 1

The idea is to bring the second order inequality to a first order inequality, using monotone operator, convexity and fixed point arguments, and will be carried out in several steps.

Let us introduce the velocity variable v = u.

The system in Problem QV is then written for a.e. t ∈ (0, T ):

                               u(t) = u 0 + t 0 v(s) ds; v(t) + A(t) v(t) + B(t) u(t) +C(t) θ (t), w -v(t) V ×V +( t 0 B(t -s) ε(u(s)) ds, ε(w) -ε(v(t))) H + ψ c (t, w) -ψ c (t, v(t)) ≥ f (t), w -v(t) V ×V ∀w ∈ V ; θ (t), η E ×E + K(t) θ (t), η E ×E + ψ thermal (t) θ (t), η E ×E = R(t)u(t), η E ×E + Q(t), η E ×E , ∀ η ∈ E; v(0) = v 0 , θ (0) = θ 0 , with the regularity v ∈ L 2 (0, T ;V ) ∩W 1,2 (0, T ;V ) ∩C(0, T ; H) θ ∈ L 2 (0, T ; E) ∩W 1,2 (0, T ; E ) ∩C(0, T ; F).
To continue, we assume in the sequel that the conditions (1.13)-(1.17) of the Theorem 1 are satisfied. Let define

W := L 2 (0, T ; H ).
We begin by Lemma 1. For all η ∈ W , there exists an unique

v η ∈ L 2 (0, T ;V ) ∩W 1,2 (0, T ;V ) ∩C(0, T ; H) satisfying              vη (t) + A(t) v η (t), w -v η (t) V ×V + (η(t), ε(w) -ε(v η (t))) H +ψ c (t, w) -ψ c (t, v η (t)) ≥ f (t), w -v η (t) V ×V , ∀ w ∈ V, a.e. t ∈ (0, T ); v η (0) = v 0 . (1.38) Moreover, ∃c > 0 such that ∀η 1 , η 2 ∈ W : v η 2 (t) -v η 1 (t) 2 H + t 0 v η 1 -v η 2 2 V ≤ c t 0 η 1 -η 2 2 H , ∀t ∈ [0, T ]. ( 1
(η = η 1 , w = v η 2 (t)), then (η = η 2 , w = v η 1 (t)).
Adding the two inequalities, we deduce that for a.e. t ∈ (0; T ):

vη 2 (t) -vη 1 (t), v η 2 (t) -v η 1 (t) V ×V + A(t) v η 2 (t) -A(t) v η 1 (t), v η 2 (t) -v η 1 (t) V ×V ≤ -(η 2 (t) -η 1 (t), ε(v η 2 (t)) -ε(v η 1 (t))) H .
Then integrating over (0,t), from (1.13)(iii) and from the initial condition on the velocity, we obtain:

∀t ∈ [0, T ], v η 2 (t) -v η 1 (t) 2 H + m A t 0 v η 2 (s) -v η 1 (s) 2 V ds ≤ - t 0 (η 2 (s) -η 1 (s), ε(v η 2 (s)) -ε(v η 1 (s))) H ds + m A t 0 v η 2 (s) -v η 1 (s) 2 H ds.
We conclude that ∃c > 0 such that ∀η 1 , η 2 ∈ W , ∀t ∈ [0, T ]:

v η 2 (t) -v η 1 (t) 2 H + t 0 v η 1 (s) -v η 2 (s) 2 V ds ≤ c t 0 η 1 (s) -η 2 (s) 2 H ds + c t 0 v η 2 (s) -v η 1 (s) 2 H ds. (1.40) Now let fix τ ∈ [0, T ]. We have ∀t ∈ [0, τ]: v η 2 (t) -v η 1 (t) 2 H ≤ c τ 0 η 1 (s) -η 2 (s) 2 H + c t 0 v η 2 (s) -v η 1 (s) 2 H ds.
Using then Gronwall's inequality, we obtain ∀τ ∈ [0, T ]:

v η 2 (τ) -v η 1 (τ) 2 H ≤ c τ 0 η 1 (s) -η 2 (s) 2 H e cT .
Finally, integrating the last inequality and reporting the result in (1.40), we get (1.39).

Here and below, we denote by c > 0 a generic constant, which value may change from lines to lines. Lemma 2. For all η ∈ W , there exists an unique

θ η ∈ L 2 (0, T ; E) ∩W 1,2 (0, T ; E ) ∩C(0, T ; F) satisfying      θη (t) + K(t) θ η (t) + ψ thermal (t) θ η (t) = R(t)v η (t) + Q(t), ∀ η ∈ E, a.e. t ∈ (0, T ); θ η (0) = θ 0 .
(1.41) Moreover, ∃c > 0 such that ∀η 1 , η 2 ∈ W :

θ η 1 (t) -θ η 2 (t) 2 F ≤ c t 0 v η 1 -v η 2 2 V , ∀t ∈ [0, T ].
(1.42)

Proof. The existence and uniqueness result verifying (1.41) follows from standard result on first order evolution equation (see e.g. [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]). Indeed we verify that from the expression of the operator R, we have Then integrating the last property over (0,t), using the strong monotonicity of K(t) and the Lipschitz continuity of R(t) : V -→ E , we deduce (1.42).

v η ∈ L 2 (0, T ;V ) =⇒ R v η ∈ L 2 (0, T ; E ), as Q ∈ L 2 (0, T ; E ) then R v η + Q ∈ L 2 (
Proof of Theorem 1.

We have now all the ingredients to prove the Theorem 1.

Consider the operator Λ : W → W defined by for all η ∈ W :

Λ η (t) = G (ε(u η (t))) + t 0 B(ts) ε(u η (s)) ds +C e (t, θ η (t)), ∀t ∈ [0, T ],

where u η (t) = u 0 + t 0 v η (s) ds, ∀t ∈ [0, T ]; u η ∈ W 1,2 (0, T ;V )∩W 2,2 (0, T ;V )∩C 1 (0, T ; H).

Then from (1.14), (1.15), and Lemma 2, we deduce that for all η 1 , η 2 ∈ W , for all t ∈ [0, T ]:

Λ η 1 (t) -Λ η 2 (t) 2 H ≤ c θ η 1 (t) -θ η 2 (t) 2 F + c t 0 v η 1 (s) -v η 2 (s) 2 V ds ≤ c t 0 v η 1 (s) -v η 2 (s) 2 V ds.
(1.43) Now using (1.43), after some algebraic manipulations, we have for any β > 0:

T 0 e -β τ Λ η 1 (τ) -Λ η 2 (τ) 2 H ≤ c β T 0 e -β τ η 1 (τ) -η 2 (τ) 2 H dτ.
We conclude from the last inequality by contracting principle that the operator Λ has a unique fixed point η * ∈ W . We verify then that the functions

  , a.e. t ∈ (0, T ), a.e. in Ω(1.15)

u(t) := u 0 + t 0 v

 0 η * , ∀t ∈ [0, T ], θ := θ η * are solutions to problem QV with the regularity (1.35), the uniqueness follows from the uniqueness in Lemma 1 and Lemma 2.

  .39) Proof. Let η ∈ W . Using[START_REF] Zeidler | Nonlinear Functional Analysis and its Applications[END_REF] II/B p. 893, we deduce the existence and uniqueness of v η . Now let η 1 , η 2 ∈ W . In (1.38) we take

  0, T ; E ). Using the assumptions(1.19) and (1.22), the operatorK(t) + ψ thermal (t) : E -→ Efor a.e. t ∈ (0, T ) is strongly monotone. Now for η 1 , η 2 ∈ W , we have for a.e. t ∈ (0; T ):θη 1 (t) -θη 2 (t), θ η 1 (t)θ η 2 (t) E ×E + K(t) θ η 1 (t) -K(t) θ η 2 (t), θ η 1 (t)θ η 2 (t) E ×E ≤ R(t) v η 1 (t) -R(t) v η 2 (t), θ η 1 (t)θ η 2 (t) E ×E .

A class of dynamic unilateral contact problems with sub-differential friction law