Alexandra Christophe-Argenvillier
email: alexandra.christophe-argenvillier@inria.fr

Stéphane Descombes
email: stephane.descombes@unice.fr

Stéphane Lanteri
email: stephane.lanteri@inria.fr

An implicit hybridized discontinuous Galerkin method for the 3D time-domain Maxwell equations

Keywords: Maxwell's equations, time-domain, implicit time stepping, discontinuous Galerkin, hybridized discontinuous Galerkin

Introduction

Generalities about the DGTD method

During the last ten years, the DGTD method has progressively emerged as a viable alternative to well established FDTD (Finite Difference Time-Domain) [START_REF] Taflove | Computational electrodynamics: the finite-difference time-domain method -3rd[END_REF] and FETD (Finite Element Time-Domain) [START_REF] Pernet | High spatial order finite element method to solve Maxwell's equations in time-domain[END_REF] methods for the numerical simulation of electromagnetic wave propagation problems in the time-domain.

The DGTD method can be considered as a finite element method where the continuity constraint at an element interface is released. While it keeps almost all the advantages of the finite element method (large spectrum of applications, complex geometries, etc.), the DGTD method has other nice properties which explain the renewed interest it gains in various domains in scientific computing:

-It is naturally adapted to a high order approximation of the unknown field.

Moreover, one may increase the degree of the approximation in the whole mesh as easily as for spectral methods but, with a DGTD method, this can also be done locally i.e. at the mesh cell level. In most cases, the approximation relies on a polynomial interpolation method but the method also offers the flexibility of applying local approximation strategies that best fit to the intrinsic features of the modeled physical phenomena.

-When the discretization in space is coupled to an explicit time integration method, the DG method leads to a block diagonal mass matrix independently of the form of the local approximation (e.g the type of polynomial interpolation). This is a striking difference with classical, continuous FETD formulations. Moreover, the mass matrix is diagonal if an orthogonal basis is chosen.

-It easily handles complex meshes. The grid may be a classical conforming finite element mesh, a non-conforming one or even a hybrid mesh made of various elements (tetrahedra, prisms, hexahedra, etc.). The DGTD method has been proven to work well with highly locally refined meshes. This property makes the DGTD method more suitable to the design of a hp-adaptive solution strategy (i.e. where the characteristic mesh size h and the interpolation degree p changes locally wherever it is needed).

-It is flexible with regards to the choice of the time stepping scheme. One may combine the discontinuous Galerkin spatial discretization with any global or local explicit time integration scheme, or even implicit, provided the resulting scheme is stable.

-It is naturally adapted to parallel computing. As long as an explicit time integration scheme is used, the DGTD method is easily parallelized. Moreover, the compact nature of method is in favor of high computation to communication ratio especially when the interpolation order is increased.

As in a classical finite element framework, a discontinuous Galerkin formulation relies on a weak form of the continuous problem at hand. However, due to the discontinuity of the global approximation, this variational formulation has to be defined at the element level. Then, a degree of freedom in the design of a discontinuous Galerkin scheme stems from the approximation of the boundary integral term resulting from the application of an integration by parts to the element-wise variational form. In the spirit of finite volume methods, the approximation of this boundary integral term calls for a numerical flux function which can be based on either a centered scheme or an upwind scheme, or a blend of these two schemes.

DGTD methods for time-domain electromagnetics

In the early 2000's, DGTD methods for time-domain electromagnetics have been studied by a few groups of researchers, most of then from the applied mathematics community. One of the most significant contributions is due to Hesthaven and Warburton [START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell's equations[END_REF] in the form of a high order nodal DGTD method formulated on unstructured simplicial meshes. The proposed formulation is based on an upwind numerical flux, nodal basis expansions on a triangle (2D case) and a tetrahedron (3D case) and a Runge-Kutta time stepping scheme. In [START_REF] Kabakian | Unstructured grid-based discontinuous Galerkin method for broadband electromagnetic simulations[END_REF], Kakbian et al. describe a rather similar approach. More precisely, the authors develop a parallel, unstructured, high order DGTD method based on simple monomial polynomials for spatial discretization, an upwind numerical flux and a fourth-order Runge-Kutta scheme for time marching. The method has been implemented with hexahedral and tetrahedral meshes. A high order DGTD method based on a strong stability preserving Runge-Kutta time scheme has been studied by Chen et al. [START_REF] Chen | High-order RKDG methods for computational electromagnetics[END_REF]. The authors also present post-processing techniques that can double the convergence order. A locally divergence-free DGTD method is formulated and studied by Cockburn et al. in [START_REF] Cockburn | Locally divergence-free discontinuous Galerkin methods for the Maxwell equations[END_REF]. In the same period, a high order nodal DGTD method formulated on unstructured simplicial meshes has also been proposed by Fezoui et al. [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF]. However, contrary to the DGTD methods discussed in [START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell's equations[END_REF] and [START_REF] Kabakian | Unstructured grid-based discontinuous Galerkin method for broadband electromagnetic simulations[END_REF], the DGTD method in [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] is non-dissipative thanks to a combination of a centered numerical flux with a second-order leap-frog time stepping scheme. The DGTD method has then been progressively considered and extended to increasinlgy more complex modeling situations by groups of researchers in the applied electromagnetics and electrical engineering communities for a wide variety of applications related to aeronautics, defense, semiconductor device fabrication, etc.

[8]-[9]-[10]- [START_REF] Alvarez | A spurious-free discontinuous Galerkin time-domain method for the accurate modeling of microwave filters[END_REF]- [START_REF] Alvarez | 3-D discontinuous Galerkin timedomain method for anisotropic materials[END_REF]- [START_REF] Dosopoulos | Non-conformal and parallel discontinuous Galerkin time domain method for Maxwell's equations: EM analysis of IC packages[END_REF]- [START_REF] Li | Cosimulation of electromagnetics-circuit systems exploiting DGTD and MNA[END_REF] to cite a few. More recently, the method has also been adopted and further developed by researchers in the nano-optics domain [START_REF] Diehl | Comparison of low-storage Runge-Kutta schemes for discontinuous Galerkin time-domain simulations of Maxwell's equations[END_REF]- [START_REF] Hille | Improving nano-optical simulations through curved elements implemented within the discontinuous Galerkin method computational[END_REF]- [START_REF] König | The Discontinuous Galerkin Time-Domain method for Maxwell's equations with anisotropic materials[END_REF]- [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF]. A full review of the nowadays numerous applications of DGDT methods would certainly require than a simple paragraph. Also worth to note, the DGTD method has been implement in commercial software such HFSS-TD (the time-domain version of the well-known HFSS software used for antenna design) [START_REF] Songoro | Keeping time with Maxwell's equations[END_REF].

Explicit versus implicit DGTD methods

From the above discussion, it is clear that the DGTD method is nowadays a very popular numerical method in the computational electromagnetics community. The works mentioned so far are mostly concerned with time explicit DGTD methods relying on the use of a single global time step computed so as to ensure stability of the simulation. It is however well known that when combined with an explicit time integration method and in the presence of an unstructured locally refine mesh, a high order DGTD method suffers from a severe time step size restriction. A possible alternative to overcome this limitation is to use smaller time steps, given by a local stability criterion, precisely where the smallest elements are located. The local character of a DG formulation is a very attractive feature for the development of explicit local time stepping schemes [START_REF] Piperno | Symplectic local time stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF]- [START_REF] Grote | Explicit local time-stepping methods for Maxwel's equations[END_REF]- [START_REF] Taube | A high-order discontinuous Galerkin method with time-accurate local time stepping for the Maxwell equations[END_REF]. An alternative approach that has been considered in [START_REF] Moya | Temporal convergence of a locally implicit discontinuous Galerkin method for Maxwell's equations[END_REF]- [START_REF] Moya | Locally implicit time integration strategies in a discontinuous Galerkin method for Maxwell's equations[END_REF] is to use a hybrid explicit-implicit (or locally implicit) time integration strategy. Such a strategy relies on a component splitting deduced from a partitioning of the mesh cells in two sets respectively gathering coarse and fine elements. In these works, a second-order explicit leap-frog scheme is combined with a second-order implicit Crank-Nicolson scheme in the framework of a non-dissipative (centered flux based) DG discretization in space. At each time step, a large linear system must be solved whose structure is partly diagonal (for those rows of the system associated to the explicit unknowns) and partly sparse (for those rows of the system associated to the implicit unknowns). The computational efficiency of this locally implicit DGTD method depends on the size of the set of fine elements that directly inluences the size of the sparse part of the matrix system. Therefore, an approach for reducing the size of the subsystem of globally coupled (i.e. implicit) unknowns is worth considering if one wants to solver very large-scale problems.

A particularly appealing solution in this context is given by the concept of hybridizable discontinuous Galerkin (HDG) method. The HDG method has been first introduced by Cockbrun et al. in [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF] for a model elliptic probelm and has been subsequently developed for a variety of PDE systems in continuum mechanics [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for partial differential equations in continuum mechanics[END_REF]. The essential ingredients of a HDG method are a local Galerkin projection of the underlying system of PDEs at the element level onto spaces of polynomials to parametrize the numerical solution in terms of the numerical trace; a judicious choice of the numerical flux to provide stability and consistency; and a global jump condition that enforces the continuity of the numerical flux to arrive at a global weak formulation in terms of the numerical trace. The HDG methods are fully implicit, high-order accurate and endowed with several unique features which distinguish themselves from other discontinuous Galerkin methods. mots importantly, they reduce the globally coupled unknowns to the approximate trace of the solution on element boundaries, thereby leading to a significant reduction in the degrees of freedom. HDG methods for the system of time-harmonic Maxwell equations have been introduced in [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF]]- [START_REF] Li | Numerical investigation of a high order hybridizable discontinuous Galerkin method for 2d time-harmonic Maxwell's equations[END_REF]- [START_REF] Li | A hybridizable discontinuous Galerkin method combined to a Schwarz algorithm for the solution of 3d time-harmonic Maxwell's equations[END_REF].

Objectives of this work

In view of devising a hybrid explicit-implicit HDG method, a preliminary step which is the focus of this work is to develop a fully implicit HDG formulation. In this paper we present such a time-implicit hybridizable discontinuous Galerkin (HDG) method for numerically solving the system of three-dimensional (3D) time-domain Maxwell equations. The method is based on a second-order Crank-Nicolson scheme for time integration whereas the discretization in space and in particular the definition of HDG traces are mostly inspired from the HDG formulations already knonw for the time-harmonic Maxwell equations [START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF]- [START_REF] Li | Numerical investigation of a high order hybridizable discontinuous Galerkin method for 2d time-harmonic Maxwell's equations[END_REF]- [START_REF] Li | A hybridizable discontinuous Galerkin method combined to a Schwarz algorithm for the solution of 3d time-harmonic Maxwell's equations[END_REF]. We study the stability of this time-implicit HDGTD (Hybridizable Discontinuous Galerkin Time-Domain) method as well as the solvability of the reduced system for the hybrid variable. The sequel of the paper is oragnized as follows: in section 2 we state the initial and boundary value problem that we want to solve, and then introduce some notations for the sequel; section 3 is the core of the paper where we describe and study the properties of the time-implicit HDGTD method; in section 4 we discuss about some implementation aspects of the method; preliminary results are presented in section 5; finally, section 6 concludes the paper. We note that the numerical results presented here before all aim at validating the proposed HDGTD method. A detailed assessment of the performances of the method for more realistic 3D problems simulated on parallel computing systems will be the subject of a subsequent paper.

Problem statement and notations

The system of 3D Maxwell equations

We consider the system of 3D time domain Maxwell's equations on a bounded polyhedral domain

Ω ⊂ R 3 ε∂ t E -curl H = 0, dans Ω × [0, T ], µ∂ t H + curl E = 0, dans Ω × [0, T ], (1)
where the symbol ∂ t denotes a time derivate, T a final time, E(x, t) and H(x, t) are the electric and magnetic fields. The dielectric permittivity tensor ε and the magnetic permeability tensor µ are varying in space, time-invariant and both positive functions. The boundary of Ω can be decomposed in ∂Ω = Γ m ∪ Γ a with Γ m ∩ Γ a = ∅ and the boundary conditions are choosen as

n × E = 0, sur Γ m × [0, T ], n × E + n × (n × H) = n × E inc + n × (n × H inc ) = g inc , sur Γ a × [0, T ]. (2)
Here n denotes the unit outward normal to ∂Ω and (E inc , H inc ) the incident field.

The first boundary condition indicates a metallic boundary condition (referring to a perfectly conducting surface), the second relation is called absorbing and takes here the form of the Silver-Müller condition. Finally, the system is supplemented with inital conditions: E 0 (x) = E(x, 0) and H 0 (x) = H(x, 0).

Notations and approximation spaces

We consider a partition T h of Ω into a set of tetraedra. Each non-empty intersection of two elements K e and K g , with e, g = {1, .., |T h |}, is called an interface. We denote by F I h the union of all interior interfaces of T h , by F B h the union of all boundary interfaces of T h , and

F h = F I h ∪ F B h .
Note that ∂T h represents all the interfaces ∂K for all K ∈ T h . As a result, an interior interface shared by two elements appears twice in ∂T h , unlike in F h where this interface is evaluated once.

For an interface F ∈ F I h , F = K e ∩ K g , let v ± be the traces of v on F from the interior of K e,g . On this interior face, we define mean values {•} and jumps • as

   {v} F = 1 2 (v + + v -), v F = n + × v + + n -× v -,
where the unit outward normal vector to K is denoted by n ± . For the boundary faces these expressions are modified as

{v} F = v + , v F = n + × v + ,
since we assume v is single-valued on the boundaries. In the following, we introduce the discontinuous finite element spaces and some basic operations on these spaces for later use. Let P pe (K e ) denote the space of polynomial functions of degree at most p e on the element K e . The discontinuous finite element space is introduced as

V h = v ∈ L 2 (Ω) 3 | v Ke ∈ [P pe (K e )] 3 , ∀K e ∈ T h .
where L 2 (Ω) is the space of square integrable functions on the domain Ω. The functions in V h are continuous inside each element but be discontinuous across the interfaces between elements. In addition, we introduce a traced finite element space M h given by

M h = η ∈ L 2 (F h ) 3 | η F f ∈ P p f (F f ) 3 , (η • n) F f = 0, ∀F f ∈ F h .
For two vectorial functions u and

v in [L 2 (D)] 3 , we denote (u, v) D = D u • v dx provided D is a domain in R 3 , and we denote < u, v > F = F u • v ds if F is a two-dimensional face.
Accordingly, for the mesh T h we have

(•, •) T h = K∈T h (•, •) K , •, • ∂T h = K∈T h •, • ∂K , •, • F h = F ∈F h •, • F , •, • Γa = F ∈F h ∩Γa •, • F . We set v t = -n × (n × v) , v n = n (n • v)
where v t and v n are the tangential and normal components of v such as v = v t + v n .

Principles and formulation of the global HDG problem

Following the classical DG approach, approximate solutions (E h ,H h ), for all t ∈ [0, T ], are seeked in the space V h × V h satisfying for all K in T h

(ε∂ t E h , v) K -(curl H h , v) K = 0, ∀v ∈ V h , (µ∂ t H h , v) K + (curl E h , v) K = 0, ∀v ∈ V h . (3)
Applying Green's formula, on both equations of (3) introduces boundary terms which are replaced by numerical traces E h and H h in order to ensure the connection between element-wise solutions and global consistency of the discretization. This leads to the formulation for all t ∈ [0, T ]

     (ε∂ t E h , v) K -(H h , curl v) K + H h , n × v ∂K = 0, ∀v ∈ V h , (µ∂ t H h , v) K + (E h , curl v) K -E h , n × v ∂K = 0, ∀v ∈ V h . (4)
It is straightforward to verify that n×v = n×v t and < H, n×v >= -< n×H, v >. Therefore, using numerical traces defined in terms of the tangential components H t h and E t h , we can rewrite (4) as

     (ε∂ t E h , v) K -(H h , curl v) K + H t h , n × v ∂K = 0, ∀v ∈ V h , (µ∂ t H h , v) K + (E h , curl v) K -E t h , n × v ∂K = 0, ∀v ∈ V h . (5)
The hybrid variable Λ h introduced in the setting of a HDG method is here defined for all the interfaces of F h as

Λ h := H t h , ∀F ∈ F h . (6)
We want to determine the fields H t h and E t h in each element K of T h by solving system [START_REF] Chen | High-order RKDG methods for computational electromagnetics[END_REF] and assuming that Λ h is known on all the faces of an element K. Let us consider a numerical trace E t h for all K given

E t h = E t h + τ K n × (Λ h -H t h ) on ∂K, (7)
where τ K is a local stabilization parameter.

Remark 1. In a classical DG method the traces of the local fields E h and H h between neighboring elements are defined as

E h = {E h } + α H H h and H h = {H h } + α E E h ,
where α H and α E are positive penalty parameters.

Remark 2. Following the HDG approach, when the hybrid variable Λ h is known for all the faces of the element K, the electromagnetic field into K e can be determined by the local system (5) using ( 6) and [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF].

Adding the contributions of ( 5) over all the elements and enforcing the continuity of the tangential component of E h , we can formulate a problem which is to find

(E h , H h , Λ h ) ∈ V h × V h × M h such that for all t ∈ [0, T ] (ε∂ t E h , v) T h -(H h , curl v) T h + Λ h , n × v ∂T h = 0, ∀v ∈ V h , (µ∂ t H h , v) T h + (E h , curl v) T h -E t h , n × v ∂T h = 0, ∀v ∈ V h , E h , η F h -Λ h , η Γa = g inc , η Γa , ∀η ∈ M h , (8)
where the last equation is called the conservativity condition with which we ask the tangential component of E h to be weakly continuous across any interface between two neighboring elements. With the definition [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF], we can apply a Green formula on the second equation of ( 8) and get for all t ∈ [0, T ]

(ε∂ t E h , v) T h -(H h , curl v) T h + Λ h , n × v ∂T h = 0, ∀v ∈ V h , (µ∂ t H h , v) T h + (curl E h , v) T h + τ n × (H h -Λ h ) , n × v ∂T h = 0, ∀v ∈ V h , n × E h , η ∂T h + τ (H t h -Λ h ) , η ∂T h -Λ h , η Γa = g inc , η Γa , ∀η ∈ M h . (9)
The main principles of the HDG method can be summarized as:

1. The DoFs of the hybrid variable are determined by solving a global linear system (from the discretization of the conservation condition) supported by the interfaces of F h ; 2. The DoFs of the electromagnetic field in each element are evaluated by solving local linear systems, more exactly for the DoFs of (E h , H h ) in the considered element.

Stability and conservation properties when

Γ a = ∅ Selecting v = E h (t)
in the first relation of ( 9) and v = H h (t) in the second relation of ( 9) and summing the two equations we obtain

1 2 d dt (ε||E h (t)|| 2 + µ||H h (t)|| 2 ) = -Λ h , n × E h (t) ∂T h -τ n × (H h (t) -Λ h ) , n × H h (t) ∂T h .
When Γ a = ∅ it follows from the third relation of ( 8) with η = Λ h that 1 2

d dt ε||E h (t)|| 2 + µ||H h (t)|| 2 = -τ n × (Λ h -H h (t)) , n × (Λ h -H h (t)) ∂T h ≤ 0, since τ > 0. Thus, the energy function E h (t) = 1 2 (ε||E h (t)|| 2 + µ||H h (t)|| 2 )
decreases in time such that E h (t) ≤ E h (0), for all t > 0. This result shows the stability in L 2 sense. In particular, this method is dissipative for the considered numerical trace for E t h in (7).

Time integration

The system of equations ( 9) can be rewriten in the form of a differential algebraic equation (DAE) such as

F (E h (t), H h (t), Λ h (t)) = 0, ( 10
)
where F is defined for all v in V h and for all η in M h by

  (ε∂ t E h , v) T h -(H h , curl v) T h + Λ h , n × v ∂T h (µ∂ t H h , v) T h + (curl E h , v) T h + τ n × (H h -Λ h ) , n × v ∂T h n × E h , η ∂T h + τ (H t h -Λ h ) , η ∂T h -Λ h , η Γa -g inc , η Γa   .
As explain in [START_REF] Hairer | The numerical solution of differential algebraic systems by Runge-Kutta methods[END_REF] the defined system is considered as a semi-explicit DAE. In this sense, the third equation of ( 10) is called algebraic equation and is considered as a constraint on the global system. The idea here is to transform the DAE into an ODE which can be straightforward solved using numerical methods. From the previews DAE, a perturbation problem is considered such as for all t ∈ [0, T ]

F E h (t), H h (t), Λ h (t) = δ Λ h (t) , (11)
with 0 < γ << 1 and

δ Λ h (t) =   0 0 γ∂ t Λ h   .
As a first step, the perturbation problem is considered

(ε∂ t E h , v) T h -(H h , curl v) T h + Λ h , n × v ∂T h = 0, ( 12
) (µ∂ t H h , v) T h + (curl E h , v) T h + τ n × (H h -Λ h ) , n × v ∂T h = 0, ( 13
)
γ∂ t Λ h + n × E h , η ∂T h + τ (H t h -Λ h ) , η ∂T h -Λ h , η Γa = g inc , η Γa . ( 14
)
From here we can discretize the problem in time. We consider a sequence of time steps

0 = t 0 < t 1 < t 2 < • • • < t N = T.
For 0 ≤ n ≤ N -1, we denote by (E n+1 h , H n+1 h , Λ n+1 h ) the numerical approximations to (E h (t n+1 ), H h (t n+1 ), Λ h (t n+1 )) at time t n+1 = (n + 1)∆t. Using a Crank-Nicolson scheme on each equation of ( 12)-( 14), the obtained system is written as:

find (E n+1 h , H n+1 h , Λ n+1 h ) ∈ V h × V h × M h such as 2ε ∆t E n+1 h , v T h -H n+1 h , curl v T h = -Λ n+1 h , n × v ∂T h + 2ε ∆t (E n , v) T h + (H n h , curl v) T h -Λ n h , n × v ∂T h , 2µ ∆t H n+1 h , v T h + curl E n+1 h , v T h = -τ n × H n+1 h -Λ n+1 h , n × v ∂T h + 2µ ∆t (H n h , v) T h -(curl E n h , v) T h -τ n × (H n h -Λ n h ) , n × v ∂T h , 2γ ∆t Λ n+1 h -Λ n h + n × E n+1 h , η ∂T h = τ (H t h ) n+1 -Λ n+1 h , η ∂T h + Λ n+1 h , η Γa + g inc,n+1 , η Γa -τ (H t h ) n -Λ n h , η ∂T h -n × E n h , η ∂T h + Λ n h , η Γa + g inc,n , η Γa .
In order to simplify the equations, we introduce ε = 2ε/∆t, μ = 2µ/∆t and

b E = 2ε ∆t (E n h , v) T h + (H n h , curl v) T h -Λ n h , n × v ∂T h , b H = 2µ ∆t (H n h , v) T h -(curl E n h , v) T h -τ n × (H n h -Λ n h ) , n × v ∂T h , b Λ = -n × E n h , η ∂T h -τ (H t h ) n -Λ n h , η ∂T h + Λ n h + g inc,n + g inc,n+1
, η Γa .

The system to be solved can be written as

εE n+1 h , v T h -H n+1 h , curl v T h + Λ n+1 h , n × v ∂T h = b E ,( 15
)
μH n+1 h , v T h + curl E n+1 h , v T h + τ n × H n+1 h -Λ n+1 h , n × v ∂T h = b H ,( 16
)
2γ ∆t Λ n+1 h -Λ n h + n × E n+1 h , η ∂T h + τ (H t h ) n+1 -Λ n+1 h , η ∂T h -Λ n+1 h , η Γa = b Λ .( 17
)
By the DAE theory γ tends to zero such as we can define the limit

lim γ→0 2γ ∆t Λ n+1 h -Λ n h = 0.
Then, the system of equations is finally given by (in order to simplify the equations, we omit the superscript n+1 )

(εE h , v) T h -(H h , curl v) T h + Λ h , n × v ∂T h = b E , (18)
(μH h , v) T h + (curl E h , v) T h + τ n × (H h -Λ h ) , n × v ∂T h = b H , (19)
n × E h , η ∂T h + τ (H t h ) -Λ h , η ∂T h -Λ h , η Γa = b Λ . (20)

Well-posedness of the local solver

For α ∈ M h , (E α h , H α h ) denotes the approximate solution at time n + 1 whose restriction to an element K ∈ T h is the solution to the local problem

(εE α h , v) K -(H α h , curl v) K + α, n × v ∂K = b E , ∀v ∈ V h (21) (μH α h , v) K + (curl E α h , v) K + τ n × (H α h -α) , n × v ∂K = b H , ∀v ∈ V h .( 22
)
Setting 21) and v = H α h in [START_REF] Taube | A high-order discontinuous Galerkin method with time-accurate local time stepping for the Maxwell equations[END_REF], from the second equation, we have

b E = b H = 0, v = E α h in (
(curl E α h , H α h ) K = -τ n × (H α h -α) , n × H α h ∂K -(μH α h , H α h ) K ,
which is inserted in the first equation and obtain

(εE α h , E α h ) K +(μH α h , H α h ) K + τ n × (H α h -α) , n × H α h ∂K + α, n × E α h ∂K = 0. ( 23
)
For α = 0 and assuming τ > 0, [START_REF] Moya | Temporal convergence of a locally implicit discontinuous Galerkin method for Maxwell's equations[END_REF] implies that E α h = H α h = 0 on K and since ε and μ are strictly positive real numbers.

Characterization of the reduced problem

Firstly, we consider equation [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF] with

Λ h = η (E η h , H η h ) and with v = E α h we obtain (εE η h , E α h ) T h -(H η h , curl E α h ) T h + η, n × E α h ∂T h = b E . (24)
Secondly, we set v = H eta h in [START_REF] Songoro | Keeping time with Maxwell's equations[END_REF] to obtain

(μH α h , H η h ) T h + (curl E α h , H η h ) T h + τ n × (H α h -α) , n × H η h ∂T h = b H . (25)
Adding ( 24) with (25) leads to

(εE η h , E α h ) T h + (μH α h , H η h ) T h + η, n × E α h ∂T h + τ n × (H α h -α) , n × H η h ∂T h = b E + b H . (26)
The conservativity condition [START_REF] Piperno | Symplectic local time stepping in non-dissipative DGTD methods applied to wave propagation problem[END_REF] for α ∈ M h where (E α h , H α h ) denotes the approximate solution, yields

n × E α h , η ∂T h + τ n × (H α h -α) , n × η ∂T h -α, η Γa = b Λ . (27)
We replace n × E α h , η ∂T h in ( 26) by ( 27) and we deduce that α is the solution of the following reduced problem for α

A h (α, η) = L h (η), ∀η ∈ M h ( 28
)
with

A h (α, η) = (εE α h , E η h ) T h + (μH α h , H η h ) T h + τ n × (α -H α h ) , n × (η -H η h ) ∂T h + α, η Γa L h (η) = b E + b H -b Λ .
Note that the bilinear form A h is symmetric positive definite. Let us now prove the existence and uniqueness of the approximate solution.

Theorem 1. If we assume τ > 0, then the HDG method ( 18)-( 20) has a unique solution (E h , H h , Λ h ) for any time iteration n ≥ 1.

Proof. We suppose that the right-hand side of ( 28) is equal to zero

A h (α, η) = 0, ∀η ∈ M h .
Next, we choose η = α and the relation becomes

A h (α, α) = ε E α h 2 +μ H α h 2 + τ n × (α -H α h ) , n × (α -H α h ) ∂T h + α, α Γa .
A sufficient condition for having A h (α, α) = 0 is to require that each of the terms are zero. The two first terms imply that E α h = H α h = 0 for all K ∈ T h since ε > 0 and μ > 0. Moreover, α = 0 on Γ a by α, α Γa = 0. Then, from H α h = 0, we have

τ n × (α -H α h ) , n × (α -H α h ) ∂T h = τ n × α, n × α ∂T h .
Assuming τ > 0 implies that α = 0 on ∂K, for all K ∈ T h . In summary, α = 0 on F h ∩ Γ a and on T h , so α = 0 everywhere.

Energy variation and unconditional stability

The total discrete electromagnetic energy in T h at time n is given by

E n = 1 2 ε||E n h || 2 + μ||H n h || 2 . ( 29
)
Lemma 1. Using ( 15)-( 17), the total discrete electromagnetic energy defined by ( 29) is non-increasing in time i.e. E n+1 ≤ E n . Then, the totally discretized problem ( 18)-( 20) is unconditionally stable.

Proof. Choosing as test functions v = (E n+1

h 16) and summing

+ E n h )/2 belonging to V h in (15) and v = (H n+1 h + H n h )/2 belonging to V h in (
1 2 ε||E n+1 h || 2 + μ||H n+1 h || 2 + Λ n+1 h , n × E n+1 h + E n h 2 ∂T h + τ n × H n+1 h -Λ n+1 h , n × H n+1 h + H n h 2 ∂T h = 1 2 ε||E n h || 2 + μ||H n h || 2 -Λ n h , n × E n+1 h + E n h 2 ∂T h -τ n × (H n h -Λ n h ) , n × H n+1 h + H n h 2 ∂T h .
The energy variation is then given by

E n+1 -E n = -Λ n+1 h + Λ n h , n × E n+1 h + E n h 2 ∂T h -τ n × H n+1 h + H n h , n × H n+1 h + H n h 2 ∂T h + τ n × Λ n+1 h + Λ n h , n × H n+1 h + H n h 2 ∂T h . From (17) with η = Λ h = Λ n+1 h +Λ n h , E h = (E n+1 h +E n h )/2 and H h = (H n+1 h +H n h )/2, we have n × E n+1 h + E n h 2 , Λ n+1 h + Λ n h ∂T h = -τ n × H n+1 h + H n h , n × H n+1 h + H n h 2 ∂T h + Λ n+1 h + Λ n h 2 .
Finally, we found the result

E n+1 -E n = -τ n × H n+1 h + H n h , n × H n+1 h + H n h 2 ∂T h -Λ n+1 h + Λ n h 2 ≤ 0.

HDG method implementation

Let us introduce some notations and definitions. For all K e ∈ T h we have

(E h Ke , H h Ke ) = (E e , H e ) : K e × K e -→ R 3 × R 3 ,
where E e (x) = E e x (x) , E e y (x) , E e z (x)

T and H e (x) = H e x (x) , H e y (x) , H e z (x) T .

The components of electromagnetic field are expanded in element K e by a linear combination of basis functions ϕ e j (x) ⊂ P pe (K e )

E e ξ (x) = N e K j=1 E e ξ [j] ϕ e j (x) , H e ξ (x) = N e K j=1 H e ξ [j] ϕ e j (x) (ξ ∈ {x, y, z}) , (30)
where N e K = (p e + 1) (p e + 2) (p e + 3) /6 the dimension of P pe , E e ξ [j] and H e ξ [j] are the DoFs of the electromagnetic field supporting by K e . For all F f ∈ F h , we have

Λ h F f = Λ f : F f -→ R 3 Λ f : (x) -→ Λ f u (x) u f + Λ f w (x) w f , (31)
where u f and w f represent coordinate axis. Λ f u and Λ f w are decomposed by linear combinations of basis functions ψ f j (x) ⊂ P p f (F f ), such as

Λ f u (x) = N f F j=1 Λ f u [j] ψ f j (x) , Λ f w (x) = N f F j=1 Λ f w [j] ψ f j (x) , (32)
where Λ f u [j] and Λ f w [j] are the DoFs of Λ f associated to F f , and

N f F = (p f + 1) (p f + 2) /2 the dimension of P p f (F f ). Remark 3. In practice, u f = u f / u f 2 and w f = w f / w f 2 with u f = n f 2 -n f 1 and w f = n f 3 -n f 1 where n f 1 , n f
∂K l e = K e ∩ K g , l ∈ {1, • • • , |ν e |}, g ∈ ν e .
We get a local indexing of faces, element per element. Moreover, we define the function σ that maps the local to the global numerotation in F h (see Figure 1)

∀F f ∈ F I h such that F f = ∂K l e ∩ ∂K k g , σ(e, l) = σ(g, k) = f, ∀F f ∈ F B h such that F f = ∂K l e ∩ (Γ m ∪ Γ m ), σ(e, l) = f.
In order to approach the values of (E h , H h) in an element K e ∈ T h , one has to know the values of the fields in each neighboring element of K e . One major criticism of DG method is that they have too many degrees of freedom (DoFs) due to nodal duplication at the element boundary interfaces

K g K e K e K g ∂K 3 e ∂K 2 e ∂K 1 g ∂K 2 g ∂K 1 e ∂K 3 g σ(e, 1) = f F f σ(g, 3) = f K e K e σ(e, 1) = f F f σ(g, 3) = f K h K g ∂K 2
|T h | e=1 6N e K ,
where N e K is the dimension of P pe (K e ), i.e. N e K = (p e + 1)(p e + 2)(p e + 3)/6. The goal of the HDG method is to reduce the number of globally coupled DoFs. For that purpose, the HDG method introduces a hybrid variable, as a new unknown, on the faces of F h . This new unknown can be determined by solving a global linear system from the discretization of a conservation condition for the hybrid variable. More exactly, the DoFs from the hybrid variable can be defined by the resolution of a global linear system. When the DoFs of the hybrid variable are known, the electromagnetic field (E h , H h ) can be locally determined element wise thanks to the resolution of local linear systems defined in each K ∈ T h . For the HDG method, the total number of globally coupled DoFs is given by

|F h | f =1 2N f F ,
where

N f F = (p f + 1) (p f + 2) /2, is the dimension of the space P p f (F f ).
Remark 4. If we consider a uniform degree of interpolation for elements and faces, i.e. p e = p f = p, ∀K e ∈ T h , ∀F f ∈ F h , then the total numbers of globally coupled DoFs are DG method:

(p + 1(p + 2)(p + 3)|T h |, HDG method: (p + 1)(p + 2)|F h |. For |F h | ≈ 2|T h |, 2/(p + 3
) is the ratio of the number of globally coupled DoFs of HDG method to DG method.

Local linear system

Taking ( 30)- [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF] and choosing the basis functions of V h as test functions in ( 18)-( 20) we obtain that E e

x , E e y , E e z , H e y , H e x , H e C (e,l) Λ σ(e,l) u Λ σ(e,l)

w = P e , (33)
A e , C (e,l) , being matrices of dimensions 6N e K × 6N e K and 6N e K × 2N

σ(e,l) F

and P e a local right-hand side.

Global linear system

A global linear system can be obtained for the DoFs of the hybrid variable Λ h in F h . Let us start by F f ∈ F I h an interior face between two elements K e and K g . This face takes the local indices l and k on K e and K g , respectively, in particular f = σ(e, l) = σ(g, k). For all η ∈ M h we have

n × E h , η ∂K l e -τ (e,l) n × (n × H h ), η ∂K l e -τ (e,l) Λ h , η ∂K l e + n × E h , η ∂K k g -τ (g,k) n × (n × H h ), η ∂K k g -τ (g,k) Λ h , η ∂K k g = 0.
(34)

Now, we consider a boundary face F f ∈ Γ a such as F f ∈ ∂K l e ∩ Γ a . For all η ∈ M h , the conservation condition takes the form

n × E h , η ∂K l e -τ (e,l) n × (n × H h ), η ∂K l e -(1 + τ (e,l) ) Λ h , η ∂K l e = g inc , η ∂K l e . (35)
Using ( 30)- [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF]

KΛ = b Λ , ( 36
)
where K is a matrix of size

|F h | f =1 2N f F × |F h | f =1 2N f F and b Λ is a vector of length |F h | f =1 2N f F .

Numerical results

The HDGTD method described in the previous sections has been programmed in Fortran 95/2003 for a renormalized version of the time-domain Maxwell equations [START_REF] Taflove | Computational electrodynamics: the finite-difference time-domain method -3rd[END_REF]. The equations are still formulated as in (1) but with ε ≡ ε r and µ ≡ µ r where ε r and µ r are relative electromagnetic parameters without dimension. Beside, in the renormalized equations, the time variable t is given in meter while the electric and magnetic fields, E and H, are obtained in Volts per meter (V/m).

Propagation of a standing wave in a cubic PEC cavity

In order to validate and to study the numerical convergence of the proposed 3D HDG method, we consider the propagation of an eigenmode in a closed cavity (Ω is the unit square) where the walls are metallic, the frequency f = √ 3/2c 0 Hz, where c 0 is the speed of light in vacuum, and the wavelength λ = 1.1547 m. The permittivity and the permeability are equals to the constant vacuum values ε = ε 0 and µ = µ 0 i.e. ε r = µ r = 1. The exact time-domaine solution is given by where the angular frequency (or pulsation) is given by ω = 2πf (rad • s -1 ). The above relations are used to initialize the E and H field components in the cavity. The parameter τ in the HDG traces is taken equal to 1. The time step is chosen as ∆t = c CFL ∆t min where ∆t min is the global minimal time step over the whole mesh. The value of the time step is chosen such that the accuracy does not affect the accuracy in space. We measure the maximal L 2 -norm of the error for a sequence of successively refined tetrahedral meshes starting from a uniform coarse mesh. The latter is obtained by subdividing a finite difference grid of the unit cube. We plot this error as a function of the square root of the number of DoFs required for the interpolation of the electric and magnetic field, in logarithmic scale. On Figure 2a, we observe an optimal convergence order for the electromagnetic field. For a too large value of the coefficient c CFL , the discretization error in time is predominant and the order is equal to 2 independently of the interpolation degree, see Figure 2b. We proceed by comparing the proposed time-implicit HDG with a classical centered flux time-implicit DGTD method [START_REF] Moya | Locally implicit time integration strategies in a discontinuous Galerkin method for Maxwell's equations[END_REF] in terms of memory occupation and CPU time for the same problem. The coefficient c CFL is set to 1. We make use of the MUMPS sparse direct solver [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF] for the inversion of the HDG system (36) and of the global linear system characterizing the time-implicit DGDT method. Simulations are performed on a workstation equipped with an Intel Xeon E5-2630@2.60 GHz processsor. The obtained results are summarized in Table 1 for two interpolation orders. We clearly see that the HDG method outperforms the DG method both on the memory requirement and CPU time metrics, especially at higher interpolation orders.

     E x (x,

Scattering of a plane wave by a dielectric sphere

We now consider the problem of the scattering of a plane wave by a dielectric sphere in free space. The computational domain Ω 1 is artificially bounded by a sphere with radius r 1 = 1.5 m on which the Silver-Müller absorbing condition is applied. The dielectric sphere Ω 2 is assumed to have a radius of r 2 = 0.5 m and bounds a material of relative permittivity ε r,2 = 2 and relative permeability µ r,2 = 1 (non-magnetic material). The medium exterior to the dielectric sphere is assumed to be vacuum, i.e. ε r,1 = µ r,1 = 1. At initial time, the electromagnetic field is set to zero while the incident plane wave is defined by

E inc x = E inc y = H inc x = H inc z = 0 and E inc z (x, y, z, t = 0) = cos(ωy), H inc y (x, y, z, t = 0) = -cos(ωy) (37)
We consider two cases. In the first one, the frequency of the incident plane wave is equal to F = 300 MHz and the total simulation time is set to T = 3.33 10 -8 s which corresponds to a propagation of the initial wave over 10 periods. In this case the time step is set to ∆t = 311∆t min where ∆t min is the smallest time step over the whole mesh which is equal to 10 -12 s (the largest time step is equal to 1.9 10 -10 s). In the second case, we set F = 600 MHz and a final time T = 1.67 10 -8 s which also corresponds to a propagation of the initial wave over 10 periods. The CFL number is here set to 155. We present results for HDG-P k methods with k = 1, 2,. The better accuracy with increasing k is here limited by the second-order time approximation (as mentioned in the previous case) and by the approximation error on the circular geometry of the interface between the two materials as well as the use of the first order Silver-Müller absorbing condition. The contour lines of the real part of the discrete Fourier transform (DFT) of the component E z evaluated are exposed in Figures 3 andFigures 4, respectively.

To assess the performances of the HDG method, we summarize in Table 2 different values, in particular the number of nonzero elements and the memory requirement for the analysis and factorization of the matrix of system (36). These value are here independent of the frequency F on the incident plane wave. Data for a simulation performed with a fully explicit centered flux-based DG method [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] are also given. Despite the overhead introduced by the resolution of the linear system for the hybrid variable (36), the time-implicit HDG method is faster in CPU time, thanks to the possibility of using a larger time step.

Conclusion

In this work we have presented a time-implicit hybridizable discontinuous Galerkin method to solve numerically the time-domain Maxwell equations in 3D. We have provided preliminary numerical results demonstrating that the method can be com- putationally efficient for moderately large 3D problems. In order to address much larger and more realistic problems, a strategy combining explicit and implicit time schemes in the spirit of [START_REF] Moya | Temporal convergence of a locally implicit discontinuous Galerkin method for Maxwell's equations[END_REF]- [START_REF] Moya | Locally implicit time integration strategies in a discontinuous Galerkin method for Maxwell's equations[END_REF] will be one of our future works. Beside, improving the accuracy in time is also an important objective. In the context of hybrid explicitimplicit time integration strategies, an original approach inspried from domain decomposition principles has recently been presented in [START_REF] Chabassier | Fourth order energy-preserving locally implicit time discretization for linear wave equations[END_REF] for the second-order wave equation. Whether such an approach can be extended to the first-order form of the system of time-domain Maxwell equations remains an open question that would be worth considering in the HDG discretization framework. Finally, devising a scalable implementation of such an hybrid explicit-implicit HDGTD method will also be a mandatory step.

2 and n f 3

 3 represent the three nodes associated to F f . Henceforth, we denote by ν e the set of indices of neighboring elements of K e . For all K e ∈ T h (e ∈ {1, • • • , |T h |), we have |ν e | faces, denoted by ∂K l e ∈ ∂T h and defined by

1 gFigure 1 :

 11 Figure 1: Function σ and indexing (local/global) in 2D: conforming mesh (left) and non-conforming mesh (right).

zT

  and Λ σ(e,l) u , Λ σ(e,l) w T , 1 ≤ l ≤ |ν e | satisfy the following system

10 1 c CFL = 1 Figure 2 :

 112 Figure 2: Standing wave in a PEC cavity: numerical convergence study.

  , replacing the test functions in (34) by the basis functions of M h

	and using (35) we can find another relation between E e x , E e y , E e z , H e y , H e x , H e z	T and
	Λ σ(e,l) u	, Λ σ(e,l)

w T , 1 ≤ l ≤ |ν e |. With (33) we obtain the following system

Table 1 :

 1 Standing wave in a PEC cavity: comparaison between the time-implicit DGTD method and the time-implicit HDGTD method.

		P 1		P 3
		DG	HDG	DG	HDG
	Matrix order	36 864	39 168	184 320	130 560
	# non-zero entries	2 988 320 1 562 112 64 448 080 17 356 800
	Fill-in ratio (%)	0.21	0.1	0.19	0.1
	Storage LU (MB)	541	115	13 442	1 271
	Overall storage (MB)	827	188	20 383	1 985
	Analysis time (s)	0.7	0.09	12	0.7
	Factorization time (s)	46	3	8 308	143
	L 2 error	4.6810 -2 4.2210 -2	3.6110 -2	3.6010 -2

Table 2 :

 2 Scattering of a plane wave by a dielectric sphere F = 300 MHz: performance results.

		Implicit	Explicit
		HDG-P 1	HDG-P 2	DG-P 1
	CFL number	311	311	0.82
	Matrix order	547 884	1 095 696	547 848
	# non-zero entries	22 671 072 90 684 288 3 456 256
	Fill-in ratio (%)	0.007	0.007	0.001
	Analysis storage (MB)	11 561	43 168	1664
	Factorization storage (MB)	8863	32 356	1648
	Analysis time (s)	3.0	7.4	0.4
	Factorization time (s)	174.9	1 725.4	0.236
	CPU time (s)	151	1 1384	16 362


Kevin Bascol

Rémi Emonet

Elisa Fromont

Jean-Marc Odobez

Unsupervised Interpretable Pattern Discovery in Time Series Using Autoencoders

We study the use of feed-forward convolutional neural networks for the unsupervised problem of mining recurrent temporal patterns mixed in multivariate time series. Traditional convolutional autoencoders lack interpretability for two main reasons: the number of patterns corresponds to the manually-fixed number of convolution filters, and the patterns are often redundant and correlated. To recover clean patterns, we introduce different elements in the architecture, including an adaptive rectified linear unit function that improves patterns interpretability, and a group-lasso regularizer that helps automatically finding the relevant number of patterns. We illustrate the necessity of these elements on synthetic data and real data in the context of activity mining in videos.

Introduction

Unsupervised discovery of patterns in temporal data is an important data mining topic due to numerous application domains like finance, biology or video analysis. In some applications, the patterns are solely used as features for classification and thus the classification accuracy is the only criterion. This paper considers different applications where the patterns can also be used for data analysis, data understanding, and novelty or anomaly detection [START_REF] Emonet | Multi-camera open space human activity discovery for anomaly detection[END_REF][START_REF] Emonet | Temporal analysis of motif mixtures using dirichlet processes[END_REF][START_REF] Du | Migration motif: a spatial -temporal pattern mining approach for financial markets[END_REF][START_REF] Sallaberry | Sequential patterns mining and gene sequence visualization to discover novelty from microarray data[END_REF]. Not all time series are of the same nature. In this work, we consider the difficult case of multivariate time series whose observations are the result of a combination of different recurring phenomena that can overlap. Examples include traffic videos where the activity of multiple cars causes the observed sequence of images [START_REF] Emonet | Temporal analysis of motif mixtures using dirichlet processes[END_REF], or aggregate power consumption where the observed consumption is due to a mixture of appliances [START_REF] Kolter | Approximate inference in additive factorial hmms with application to energy disaggregation[END_REF]. Unlike many techniques from the data mining community, our aim is not to list all recurrent patterns in the data with their frequency but to reconstruct the entire temporal documents by means of a limited and unknown number of recurring patterns together with their occurrence times in the data. In this view, we want to un-mix multivariate time series to recover how they can be decomposed in terms of recurrent temporally-structured patterns. Following the conventions used in [START_REF] Emonet | Temporal analysis of motif mixtures using dirichlet processes[END_REF], we will call a temporal pattern a motif, and an input multivariate time series a temporal document.

Artificial neural networks (or deep learning architectures) have (re)become tremendously popular in the last decade due to their impressive, and so far not beaten, results in image classification, speech recognition and natural language processing. In particular, autoencoders are artificial neural networks used to learn a compressed, distributed representation (encoding) for a set of data, typically for the purpose of dimensionality reduction. It is thus an unsupervised learning method whose (hidden) layers contain representations of the input data sufficiently powerful for compressing (and decompressing) the data while loosing as few information as possible. Given the temporal nature of your data, our pattern discovery task is fundamentally convolutional (the same network is applied at any instant and is thus time-shift invariant) since it needs to identify motifs whatever their time(s) of occurrence. To tackle this task, we will thus focus on a particular type of autoencoders, the convolutional ones. However, while well adapted for discriminative tasks like classification [START_REF] Baccouche | Spatio-Temporal Convolutional Sparse Auto-Encoder for Sequence Classification[END_REF], the patterns captured by (convolutional) autoencoders are not fully interpretable and often correlated.

In this paper, we address the discovery of interpretable motifs using convolutional auto-encoders and make the following contributions:

we show that the interpretability of standard convolutional autoencoders is limited; we introduce an adaptive rectified linear unit (AdaReLU) which allows hidden layers to capture clear occurrences of motifs, we propose a regularization inspired by group-lasso to automatically select the number of filters in a convolutional neural net, we show, through experiments on synthetic and real data, how these elements (and others) allow to recover interpretable motifs 3 . It is important to note that some previous generative models [START_REF] Varadarajan | A sequential topic model for mining recurrent activities from long term video logs[END_REF][START_REF] Emonet | Temporal analysis of motif mixtures using dirichlet processes[END_REF] have obtained very good results on this task. However, their extensions to semi-supervised settings (i.e. with partially labelled data) or hierarchical schemes are cumbursome to achieve. In contrast, in this paper, to solve the same modeling problem we present a radically different method which will lend itself to more flexible and systematic end-to-end training frameworks and extensions.

The paper is organized as follows. In Section 2, we clarify the link between our data mining technique and previous work. Section 3 gives the details of our method while Section 4 shows experiments both on synthetic and real data. We conclude and draw future directions in Section 5.

Related Work

Our paper shows how to use a popular method (autoencoders) to tackle a task (pattern discovery in time series) that has seldom been considered for this type of method. We thus briefly review other methods used in this context and then, other works that use neural networks for unsupervised time series modeling.

Unsupervised pattern discovery in time series. Traditional unsupervised approaches that deal with time series do not aim at modeling series but rather at extracting interesting pieces of the series that can be used as high level descriptions for direct analysis or as input features for other algorithms. In this category fall all the event-based (e.g. [START_REF] Peng | Mining temporal patterns in time intervalbased data[END_REF][START_REF] Zhou | Unsupervised temporal commonality discovery[END_REF][START_REF] Marwah | A temporal motif mining approach to unsupervised energy disaggregation: Applications to residential and commercial buildings[END_REF]), sequence [START_REF] Mooney | Sequential pattern mining -approaches and algorithms[END_REF] and trajectory mining methods [START_REF] Zheng | Trajectory data mining: An overview[END_REF]. On the contrary of the previously cited methods, we do not know in advance the occurrence time, type, length or number of (possibly) overlapping patterns that can be used to describe the entire multivariate time series. These methods cannot be directly used in our application context.

The generative methods for modeling time series assume an apriori model and estimate its parameters. In the precursor work of [START_REF] Oates | Peruse: An unsupervised algorithm for finding recurring patterns in time series[END_REF], the unsupervised problem of finding patterns was decomposed into two steps, a supervised step involving an oracle who identifies patterns and series containing such patterns and an EM-step where a model of the series is generated according to those patterns. In [START_REF] Mehta | Funcica for time series pattern discovery[END_REF], the authors propose a functional independent component analysis method for finding linearly varying patterns of activation in the data. They assume the availability of pre-segmented data where the occurrence time of each possible pattern is known in advance. Authors of [START_REF] Kolter | Approximate inference in additive factorial hmms with application to energy disaggregation[END_REF] address the discovery of overlapping patterns to disaggregate the energy level of electric consumption. They propose to use additive factorial hidden Markov models, assuming that the electrical signal is univariate and that the known devices (each one represented by one HMM) have a finite known number of states. This also imposes that the motif occurrences of one particular device can not overlap. The work of [START_REF] Emonet | Temporal analysis of motif mixtures using dirichlet processes[END_REF] proposes to extract an apriori unknown number of patterns and their possibly overlapping occurrences in documents using Dirichlet processes. The model automatically finds the number of patterns, their length and occurrence times by fitting infinite mixtures of categorical distributions to the data. This approach achieved very good results, but its extensions to semi-supervised settings [START_REF] Tavenard | Time-sensitive topic models for action recognition in videos[END_REF] or hierarchical schemes [START_REF] Chockalingam | Localized anomaly detection via hierarchical integrated activity discovery[END_REF] were either not so effective [START_REF] Tavenard | Time-sensitive topic models for action recognition in videos[END_REF] or more cumbursome [START_REF] Chockalingam | Localized anomaly detection via hierarchical integrated activity discovery[END_REF]. In contrast, the neural network approach of this paper will lend itself to more flexible and systematic end-to-end training frameworks and extensions.

Networks for time series mining. A recent survey [START_REF] Karlsson | A review of unsupervised feature learning and deep learning for time-series modeling[END_REF] reviews the networkbased unsupervised feature learning methods for time series modeling. As explained in Sec. 1, autoencoders [START_REF] Ranzato | Efficient learning of sparse representations with an energy-based model[END_REF] and also Restricted Boltzmann Machines (RBM) [START_REF] Hinton | [END_REF] are neural networks designed to be trained from unsupervised data. The two types of networks can achieve similar goals but differ in the objective function and related optimization algorithms. Both methods were extended to handle time series [START_REF] Memisevic | Unsupervised learning of image transformations[END_REF][START_REF] Baccouche | Spatio-Temporal Convolutional Sparse Auto-Encoder for Sequence Classification[END_REF], but the goal was to minimize a reconstruction error without taking care of the interpretability or of finding the relevant number of patterns. In this paper, we show that convolutional autoencoders can indeed capture the spatio-temporal structure in temporal documents. We build on the above works and propose a model to discover the right number of meaningful patterns in the convolution filters, and to generate sparse activations.

Motif Mining with Convolutional Autoencoders (AE)

Convolutional AEs [START_REF] Masci | Stacked convolutional autoencoders for hierarchical feature extraction[END_REF] are particular AEs whose connection weights are constrained to be convolution kernels. In practice, this means that most of the learned parameters are shared within the network and that the weight matrices which store the convolution filters can be directly interpreted and visualized. Below, we first present the traditional AE model and then introduce our contributions to enforce at the same time a good interpretability of the convolutional filters and a clean and sparse activation of these filters.

Classical Convolutional Autoencoders

A main difference between an AE and a standard neural network is the loss function used to train the network. In an AE, the loss does not depend on labels, it is the reconstruction error between the input data and the network output. Fig. 1 illustrates the main network modeling components of our model. In our case, a training example is a multivariate time series x whose L time steps are described by a vector x (:,t) ∈ R d , and the network is parameterized by the set of weights W = { e W, d W} involved in the coding and decoding processes. If we denote by X = {x b ∈ R L×d , b = 1 . . . N } the set of all training elements, the estimation of these weights is classically conducted by optimizing the cost function C(W, X) = M SE(W, X)+R reg (W, X) where the Mean Squared Error (MSE) reconstruction loss can be written as:

M SE(W, X) = 1 N N b=1 d i=1 L t=1 x b (i,t) -o b (i,t) 2 (1)
where o b (which depends on parameters W) is the AE output of the b th input document. To avoid learning trivial and unstable mappings, a regularization term R reg is often added to the MSE and usually comprises two terms. The first one, known as weight decay as it avoids unnecessary high weight values, is a 2 norm on the matrix weights. The second one (used with binary activations) consists of a Kullback-Leibler divergence M j=1 KL(ρ||ρ j ) encouraging all hidden activation units to have their probability of activation ρj estimated across samples to be close to a chosen parameter ρ, thus enforcing some activation sparsity when ρ is small. The parameters are typically learned using a stochastic gradient descent algorithm (SGD) with momentum using an appropriate rate scheduling [START_REF] Darken | Note on learning rate schedules for stochastic optimization[END_REF].

Interpretable Pattern Discovery with Autoencoders

In our application, the learned convolution filters should not only minimize the reconstruction error but also be directly interpretable. Ideally, we would like to only extract filters which capture and represent interesting data patterns, as illustrated in Fig. 2-c-d. To achieve this, we add a number of elements in the network architecture and in our optimization cost function to constrain our network appropriately.

Enforcing non-negative decoding filters. As the AE output is somehow defined as a linear combination of the decoding filters, then these filters can represent the patterns we are looking for and we can interpret the hidden layers activations a (see Fig. 1) as the occurrences of these patterns. Thus, as our input is non-negative (a temporal document), we constraint the decoding filters weights to be non-negative by thresholding them at every SGD iteration. The assumption that the input is non-negative holds in our case and it will also hold in deeper AEs provided that we use ReLU-like activation functions. Note that for encoding, we do not constrain filters so they can have negative values to compensate for the pattern auto-correlation (see below).

Sparsifying the filters. The traditional 2 regularization allows many small but non-zero values. To force these values to zero and thus get sparser filters, we replaced the 2 norm by the sparsity-promoting norm 1 known as lasso:

R las (W) = M f =1 d i=1 L f k=1 e W f (i,k) + M f =1 d i=1 L f k=1 d W f (i,k) (2)
Encouraging sparse activations. The traditional KL divergence aims at making all hidden units equally useful on average, whereas our goal is to have the activation layer to be as sparse as possible for each given input document. We achieve this by encouraging peaky activations, i.e. of low entropy when seen as a document-level probability distribution, as was proposed in [START_REF] Varadarajan | A sparsity constraint for topic models -application to temporal activity mining[END_REF] when dealing on topic models for motif discovery. This results in an entropy-based regularization expressed on the set A = {a b } of document-level activations:

Rent(A) = - 1 N N b=1   M f =1 L-L f +1 t=1 âb f,t log âb f,t   with âb f,t = a b f,t M f =1 L-L f +1 t=1 a b f,t (3)
Local non-maximum activation removal. The previous entropy regularizer encourages peaked activations. However, as the encoding layer remains a convolutional layer, if a filter is correlated in time with itself or another filter, then the activations cannot be sparse. This phenomenon is due to the feed forward nature of the network, where activations depend on the input, not on each others: hence, no activation can inhibit its neighboring activations. To handle this issue we add a local non-maximum suppression layer which, from a network perspective, is obtained by convolving activations with a temporal Gaussian filter, subtracting from the result the activation intensities, and applying a ReLU, focusing in this way spread activations into central peaks.

Handling distant filter correlations with AdaReLU. The Gaussian layer cannot handle non local (in time) correlations. To handle this, we propose to replace the traditional ReLU activation function by a novel one called adaptive ReLU. AdaReLU works on groups of units and sets to 0 all the values that are below a percentage (e.g., 60%) of the maximal value in the group. In our architecture, AdaReLU is applied separately on each filter activation sequence.

Finding the true number of patterns. One main advantage and contribution of our AE-based method compared to methods presented in Section 2 is the possibility to discover the "true" number of patterns in the data. One solution to achieve this is to introduce in the network a large set of filters and "hope" that the learning leads to only a few non null filters capturing the interesting patterns. However, in practice, standard regularization terms and optimizations tend to produce networks "using" all or many more filters than the number of true patterns which results in partial and less interpretable patterns. To overcome this problem, we propose to use a group lasso regularization term called 2,1 norm [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] that constrains the network to "use" as few filters as possible. It can be formulated for our weight matrix as:

Rgrp(W) = M f =1 d i=1 L f k=1 e W f (i,k) 2 + M f =1 d i=1 L f k=1 d W f (i,k) 2 (4)
Overall objective function. Combining equations ( 1), ( 2), ( 4) and (3), we obtain the objective function that is optimized by our network:

C(W, X) = M SE(W, X) + λ las R las (W) + λgrpRgrp(W) + λentRent(A(W, X)) (5)
4 Experiments

Experimental Setting

Datasets. To study the behavior of our approach, we experimented with both synthetic and real video datasets. The synthetic data were obtained using a known generation process: temporal documents were produced by sampling random observations of random linear combinations of motifs along with salt-andpepper noise whose amount was defined as a percentage of the total document intensities (noise levels: 0%, 33%, 66%). Six motifs (defined as letter sequences for ease of visualization) were used. A document example is shown in Fig. 2-a, where the the feature dimension (d = 25) is represented vertically, and time horizontally (L = 300). For each experiments, 100 documents were generated using this process and used to train the autoencoders. This controlled environment allowed us to evaluate the importance of modeling elements. In particular, we are interested in i) the number of patterns discovered (defined as the non empty decoding filters 4 ; ii) the "sharpness" of the activations; and iii) the robustness of our method according to parameters like λ lasso , λ grp , λ ent , the number of filters M , and the noise level.

We also applied our approach on videos recorded from fixed cameras. We used videos from the QMUL [START_REF] Hospedales | A markov clustering topic model for mining behavior in video[END_REF] and the far-field datasets [START_REF] Varadarajan | A sequential topic model for mining recurrent activities from long term video logs[END_REF]. The data pre-processing steps from the companion code of [START_REF] Emonet | Temporal analysis of motif mixtures using dirichlet processes[END_REF] were applied. Optical flow features were obtained by estimating, quantifying, and locally collecting optical flow over 1 second periods. Then, temporal documents were obtained by reducing the dimensionality of these to d = 100, and by cutting videos into temporal documents of size L = 300 time steps.

Architecture details and parameter setting. The proposed architecture is given in Fig. 1. As stated earlier, the goal of this paper is to make the most of a convolutional AE with a single layer (corresponding to the activation layer) 5 . Weights are initialized according to a uniform distribution between 0 and 1 d * L f . In general, the filter length L f should be large enough to capture the longest expected recurring pattern of interest in the data. The filter length has been set to L f = 45 in synthetic experiments, which is beyond the longer motif of the ground-truth. In the video examples, we used L f = 11, corresponding to 10 seconds, and which allows to capture the different traffic activities and phases of our data [START_REF] Varadarajan | A sequential topic model for mining recurrent activities from long term video logs[END_REF].

Results on the Synthetic Dataset

Since we know the "true" number of patterns and their expected visualization, we first validate our approach by showing (see Fig. 2-c) that we can find a set of parameters such that our filters exactly capture our given motifs and the number of non empty filters is exactly the "true" number of motifs in the dataset even when this dataset is noisy (this is also true for a clean dataset). In this case (see Fig. 2-e) the activations for the complete document are, as expected, sparse and "peaky". The output document (see Fig. 2-b) is a good un-noisy reconstruction of the input document shown in Fig. 2-a.

In Fig. 3, we evaluate the influence of the given number of filters M and the noise level on both the number of recovered motifs an the MSE while fixing the parameters as in Fig. 2. We can see that with this set of parameters, the AE is able to recover the true number of filters for the large majority of noise levels and values of M . For all noise levels, we see from the low MSE that the AEs is able to well reconstruct the original document as long as the number of given filters is at least equal to the number of "true" patterns in the document.

Model selection: influence of λ grp . Fig. 4 shows the number of non zero filters in function of λ grp and of the noise level for the synthetic dataset with 6 known motifs when using 12 filters (left) and 16 filters (right). The light blue area is the area in which the AEs was able to discover the true number of patterns. With no group lasso regularization (λ grp = 0), the AE systematically uses all the available filters capturing the original patterns (see 2 nd , 4 th or 5 th filters in Fig. 2-d), redundant variants of the same pattern (filters 1 st and 3 rd in Fig. 2-d) or a more difficult to interpret mix of the patterns (filters 6 th and 7 th in Fig. 2-d).

On the contrary, with too high values of λ grp , the AE does not find any patterns (resulting in a high MSE). A good heuristic to set the value of λ grp could thus be to increase it as much as possible until the resulting MSE starts increasing.

In the rest of the experiments, λ grp is set equal to 2.

Influence of λ ent , λ lasso , AdaReLU, and Non-Local Maxima suppression. We have conducted the same experiments as in Fig. 2 on clean and noisy datasets (up to 66% of noise) with M =3, M =6 M =12 to assess the behavior of our system when canceling the parameters: 1) λ ent that controls the entropy of the activation layer, 2) λ las , the lasso regularizer 3) the AdaReLU function (we used a simple ReLU in the encoding layer instead) and 4) the Non-Local Maxima activation suppression layer. In all cases, all parameters but one were fixed according to the best set of values given in Fig. 2. For lack of space, we do not give all the corresponding figures but we comment the main results. The λ ent is particularly important in the presence of noise. Without noise and when this parameter is set to 0, the patterns are less sharp and smooth and the activations are more spread along time with much smaller intensities. Fig. 4: Evolution of the number of non-zero filters (sparsity) with respect to the noise level when we vary the parameter λ grp (λ glasso in the figure) that controls the group lasso regularization for the synthetic dataset with 6 known motifs when using 12 filters (right) and 16 filters (left).

However, the MSE is as low as for the default parameters. In the presence of noise (see Fig. 2-f), the AE is more likely to miss the recovery of some patterns even when the optimal number of filters is given (e.g. in some experiments only 5 out of the 6 filters were not empty) and the MSE increases a lot compared to experiments on clean data. This shows again that the MSE can be a good heuristic to tune the parameters on real data. The λ las has similar effects with and without noise: it helps removing all the small activation values resulting in much sharper (and thus interpretable) patterns. The non-local maximum suppression layer (comprising the Gaussian filter) is compulsory in our proposed architecture. Indeed, without it, the system was not able to recover any patterns when M =3 (and only one blurry "false" pattern in the presence of noise). When M =6, it only captured 4 patterns (out of 6) in the clean dataset and did not find any in the noisy ones. When M =12, it was able to recover the 6 original true patterns in the clean dataset but only one blurry "false" pattern in the noisy ones.

The AdaReLU function also plays an important role to recover interpretable patterns. Without it (using ReLU instead) the patterns recognized are not the "true" patterns, they have a very low intensity and are highly auto-correlated (as illustrated by the activations in Fig. 2-g).

Results on the Real Video Dataset

Due to space limitations, we only show in Fig. 5 some of the obtained results. The parameters were selected using grid search by minimizing the MSE on the targeted dataset. For instance, on the Junction 1 dataset, the final parameters used are λ las = 0.2, λ grp = 50, λ ent = 5. Note that this is larger than in the synthetic case but the observation size is also much larger (100 vs 25) and the filters are thus sparser in general. In the Junction 1 dataset, the autoencoder recovers 4 non-empty and meaningful filters capturing the car activities related to the different traffic signal cycles, whereas in the far-field case, the main trajectories of cars were recovered as also reported in [START_REF] Varadarajan | A sequential topic model for mining recurrent activities from long term video logs[END_REF].

Conclusion

We have shown that convolutional AEs are good candidate unsupervised data mining tools to discover interpretable patterns in time series. We have introduced a number of layers and regularization terms to the standard convolutional AEs to enforce the interpretability of both the convolutional filters and the activations in the hidden layers of the network. The filters are directly interpretable as spatio-temporal patterns while the activations give the occurrence times of each patterns in the temporal document. This allow us to un-mix multivariate time series. A direct perspective of this work is the use of multi-layer AEs to capture combination of motifs. If this was not the aim of this article, it may help to reduce the number of parameters needed to obtain truly interpretable patterns and capture more complex patterns in data.

Fig. 1 :

 1 Fig. 1: Autoencoder architecture. Temporal documents of L time steps of d dimensional observations are encoded (here using M convolutional filters of size d × L f forming the e W weights) to produce an activation layer. A decoding process (symmetric to encoding; parameterized by the weights d W of M decoding convolutional filters of size d × L f ) regenerates the data.

Fig. 2 :

 2 Fig. 2: Results on the synthetic data built from 6 patterns, with 66% of noise, M = 12 filters, and unless stated otherwise λ las = 0.004, λ grp = 2, λ ent = 0.2. a) Sample document; b) Obtained output reconstructed document; c) Weights of seven out of the 12 obtained filters (the 5 remaining filters are empty); d) Weights of seven filters when not using group lasso, i.e. with λ grp = 0 (note that the 5 remaining filters are non empty); (e,f,g) Examples of activation intensities (colors correspond to a given filter) with default parameters (e); without the entropy sparsifying term (λ ent = 0) (f); with ReLU instead of AdaReLU (g).

Fig. 3 :

 3 Fig. 3: Influence of the given number of filters M and the noise level (0%, 33% and 66%) on: a) the number of recovered motifs and b) the Mean Squared Error. Experiments on the synthetic dataset with λ las = 0.004, λ grp = 2, λ ent = 0.2.

Fig. 5 :

 5 Fig. 5: Traffic patterns. M =10 filters. a) The four motifs recovered on the Junction 1 dataset, (6 empty ones are not shown). b) Two filters (out of the five recovered) on the far-field dataset.

We consider a filter empty if the sum of its weights is lower or equal to 1 2 (the average sum value after initialization).

Note however that the method can be generalized to hierarchical motifs using more layers, but then the interpretation of results would slightly differ.

Acknowledgement

This work has been supported by the ANR project SoLStiCe (ANR-13-BS02-0002-01).