Juliette Chabassier

Sébastien Imperiale

Introduction and study of fourth order theta schemes for linear wave equations

Keywords: Wave equations, High order numerical methods, Time discretization, Thetascheme, Modified equation Équations d'ondes, méthodes numériques d'ordre élevé, discrétisation en temps, theta schéma, équation modifiée Inria

A new class of high order, implicit, three time step schemes for semi-discretized wave equations is introduced and studied. These schemes are constructed using the modified equation approach, generalizing the θ-scheme. Their stability properties are investigated via an energy analysis, which enables us to design super convergent schemes and also optimal stable schemes in terms of consistency errors. Specific numerical algorithms for the fully discrete problem are tested and discussed, showing the efficiency of our approach compared to second order θ-schemes.

Introduction et analyse de theta schémas d'ordre quatre pour les équations d'ondes linéaires.

Introduction

Linear wave equations play a great role in scientific modeling and are present in many fields of physics. For instance, they arise in the Maxwell equations, the acoustic equation and the elastodynamic equation. A discrete approximation of their solutions can be found with numerical simulations. Spatial discretization of the above equations using classical finite elements methods often leads to a semi-discretized problem of the form: find

u h ∈ C 2 (t, R N ) such that M h d 2 dt 2 u h + K h u h = 0, u h (0) = u 0,h , du h dt (0) = u 1,h , (1)
where u h (t) is a vector-unknown in R N , M h a symmetric positive definite matrix and K h a symmetric positive semi-definite matrix.

Several approaches can be adopted to tackle the time discretization of problem (1). The so called "conservative methods" (as for instance the leap frog scheme) preserve a discrete energy which is consistent with the physical energy. They can be shown to be stable as soon as some positivity properties of the discrete energy are satisfied, which generally imposes a restriction, known as the CFL condition, on the time step depending on the matrices M h and K h . The leap frog scheme enters a more general class of three points time step, energy preserving, implicit schemes called θ-schemes, which are parametrized by a real number θ. The over cost of these implicit schemes compared to explicit ones is balanced by the fact that stability conditions allow for a bigger time step.

For simple configurations with simple finite elements methods (such as P 1 triangular elements), explicit schemes show good performances. However they have two major drawbacks in complex configurations that have not yet been completely solved:

• If the mesh has different scales of elements, or if the equations involve variable coefficients with strong contracts, the time step must be adapted to the worst situation (for instance the smallest element) because of the CFL condition. A natural way to avoid this restriction is to use local time stepping techniques which divide into two categories. The locally implicit technique, as developed in [START_REF] Collino | A conservative space-time mesh refinement method for the 1-d wave equation. Part I: Construction[END_REF], [START_REF] Collino | A conservative space-time mesh refinement method for the 1-D wave equation. II. Analysis[END_REF], [START_REF] Kampanis | Effective computational methods for wave propagation[END_REF] and [START_REF] Becache | Space-time mesh refinement for elastodynamics[END_REF], is optimal in term of CFL restriction but "only" second order accurate in time, and requires the inversion of interface matrices. The fully explicit local time stepping, as developed in [START_REF] Diaz | Energy conserving explicit local time-stepping for second-order wave equations[END_REF], achieves higher order time stepping but without (up to now) a full control over the CFL condition.

• If the mass matrix is non diagonal or non block-diagonal, its inversion (at least one time per iteration) can lead to a dramatic over cost of the explicit schemes, whereas no over cost is observed with implicit schemes (see remark 4.1).

The extension of conservative time discretization schemes to higher orders of accuracy is a natural question. A popular way to design explicit high order three points schemes is the modified equation approach. In this article we extend this approach to design new high order implicit schemes which are stable and present some optimal properties.

The paper is organized as follows. Section 2 recalls some well-known results concerning the leap frog and θ-schemes. In a conclusive remark we present, in the very simple case of the θ-scheme, the non standard approach that we will choose to follow in the rest of the paper. In section 3 we construct a family of energy preserving implicit fourth order schemes, parametrized by two real numbers (θ, ϕ). In section 4 we discuss the existence of their discrete solution and some practical aspects of computation which can reduce numerical cost. Section 5 is devoted to the study of the stability of the newly introduced schemes via energy techniques. The search for "optimal"

Inria schemes is presented in section 6, it is done by adjusting (θ, ϕ) to increase accuracy. Finally, numerical results compare theses schemes to classical schemes in section 7.

In the following we will consider the semi-discretized problem

d 2 dt 2 u h + A h u h = 0, u h (0) = u 0,h , du h dt (0) = u 1,h , (2)
with A h a symmetric positive semi-definite matrix. With no loss of generality: the analysis done below is valid if A h = M -1 h K h . Extensive use will be made of the spectral radius of the matrix A h , defined as ρ(A h ) = sup v =1 A h v• v and coinciding with the greatest eigenvalue of A h . Its exact expression can be given in simple cases, as for instance the 1D wave equation with constant speed c, using finite differences on a regular mesh of size h, for which ρ(A h ) = 4c 2 /h 2 . In other cases, the cost of its numerical evaluation (for example with the power iteration method) is negligible compared to the numerical resolution of the equation.

Classical results

In this section we recall the definitions and some properties of the classical leap-frog and θschemes, which are widely used for the time discretization of wave equations. In the following we denote ∆t > 0 the time step of the numerical method.

Preliminary notations

The centered second order approximation of the second order derivative in time of any function t → f (t) will be denoted

D 2 ∆t f (t) = f (t + ∆t) -2 f (t) + f (t -∆t) ∆t 2 . ( 3
)
Assuming infinite smoothness on f , let us use a Taylor expansion to write the truncation error of the previous quantity:

D 2 ∆t f (t) = d 2 dt 2 f (t) + 2 ∞ m=1 ∆t 2m (2m + 2)! d 2m+2 dt 2m+2 f (t). (4)
Classical θ-schemes are based upon the use of a three points centered approximation of f (t) which, for θ ∈ R, is defined by

{f (t)} θ = θ f (t + ∆t) + (1 -2θ) f (t) + θ f (t -∆t). (5)
Assuming again infinite smoothness on f , we can write the truncation error of this new quantity:

{f (t)} θ = f (t) + 2θ ∞ m=1 ∆t 2m (2m)! d 2m dt 2m f (t). (6)
Both the leap frog scheme and the θ-scheme use finite differences to discretize time in order to compute an approximation of the semi discrete solution u h of (2). Consequently, the unknowns of those schemes stand for the values of u h at time t n = n ∆t :

u n h ≃ u h (t n ).
The discrete versions of (3) and (5), using the same symbols, are

D 2 ∆t u n h = u n+1 h -2u n h + u n-1 h ∆t 2 , {u n h } θ = θ u n+1 h + (1 -2θ)u n h + θ u n-1 h . (7)
The following algebraic relations will be useful: for all θ ∈ R,

{u n h } θ = θ θ ′ {u n h } θ ′ + θ ′ -θ θ ′ u n h , ∀ θ ′ = 0, (8)
{u n h } θ = u n h + θ∆t 2 D 2 ∆t u n h , (9)
{u n h } θ = (θ - 1 4 )∆t 2 D 2 ∆t u n h + {u n h } 1/4 . ( 10
)

Classical θ-schemes

The second order accurate leap-frog scheme reads

D 2 ∆t u n h + A h u n h = 0, (11)
which is stable under a restriction on the time step called CFL condition. More precisely, it is possible to show using energy techniques that stability requires the following inequality to hold:

∆t 2 ≤ 4 ρ(A h ) , (12)
where ρ(A h ) is the spectral radius of the matrix A h . This scheme enters a more general class of schemes called θ-schemes (convergence proofs are given in [START_REF] Karaa | Finite element theta-schemes for the acoustic wave equation[END_REF] and some alternatives are studied in [START_REF] Belytschko | Stability of explicit-implicit mesh partitions in time integration[END_REF], [START_REF] Rylander | Stability of Explicit-Implicit Hybrid Time-Stepping Schemes for Maxwell's Equations[END_REF] and [START_REF] Hui Liang | Stability of theta-schemes in the numerical solution of a partial differential equation with piecewise continuous arguments[END_REF]) :

D 2 ∆t u n h + A h {u n h } θ = 0. ( 13
)
The leap-frog scheme corresponds to the choice θ = 0. Other choices lead to implicit schemes. When θ < 1/4, these schemes are stable under the CFL condition

∆t 2 ≤ 4 (1 -4θ)ρ(A h ) . ( 14
)
As long as θ ≥ 1/4 they are unconditionally stable.

Let u h (t) be solution of (2). The truncation error of schemes (13) up to fourth order can be written

D 2 ∆t u h (t n ) + A h {u h (t n )} θ = ∆t 2 1 12 d 2 dt 2 + θA h d 2 dt 2 u h (t n ) + O(∆t 4 ). (15)
Knowing that u h (t) is solution of (2), we can replace the second order time derivative by the matrix -A h , giving

D 2 ∆t u h (t n ) + A h {u h (t n )} θ = -∆t 2 θ - 1 12 A 2 h u h (t n ) + O(∆t 4 ). (16)
This expression shows that we obtain second order accuracy except for θ = 1/12, which cancels out the first consistency error term, giving a fourth order scheme.

In some situations, the time step ∆t is related to the physics (for sampling reasons) and the Inria spectral radius ρ(A h ) (via the spatial discretization) is related to the geometry or the a priori knowledge of the solution (in the case of rapidly varying coefficients or boundary layers). Therefore, they have to be considered as data for the numerical analyst. It is then reasonable to choose a numerical scheme suited to these parameters instead of choosing discretization parameters suited to the problem. A good criterion to compare numerical schemes is the coefficient of the first term of the consistency error, which is given, for the specific case of the θ-schemes, by (θ-1/12). As these schemes give a degree of freedom, adapting θ to the product ∆t 2 ρ(A h ) makes it possible to decrease the consistency error while providing a stable scheme. This "optimization" problem is straightforward here: two cases arise, either ∆t 2 ρ(A h ) ≤ 6, in which case the value θ = 1/12 leads to a fourth order stable scheme (relation (14) holds), or if ∆t 2 ρ(A h ) > 6, the choice of θ that provides a stable scheme and minimizes the consistency error is

θ = 1 4 - 1 ∆t 2 ρ(A h ) . ( 17
)
The non standard approach mentioned above will be adopted in the discussion that follows: instead of choosing a time step ∆t that provides a good accuracy on a given mesh with a given numerical scheme, we will invert the reasoning and choose the best numerical scheme for a given couple (∆t, ρ(A h )).

3 Construction of a family of fourth order implicit schemes

Modified equation technique for the leap frog scheme

The "modified equation" technique, introduced in [SB87], enables to construct higher order schemes from the second order leap frog scheme. More precisely, the order 2p scheme is obtained by adding to the leap frog schemes terms that compensate the (p -1) first terms of the truncation error. This error reads, for u h (t) solution of (2):

D 2 ∆t u h (t n ) + A h u h (t n ) = 2 ∞ m=1 ∆t 2m (2m + 2)! d 2m+2 dt 2m+2 u h (t n ). (18)
Again, we can replace the second order time derivative by the matrix -A h :

D 2 ∆t u h (t n ) + A h u h (t n ) = 2 ∞ m=1 (-1) p+1 ∆t 2m (2m + 2)! A m+1 h u h (t n ). (19)
Truncating the series up to the first term, and approaching u h (t n ) with u n h gives the fourth order scheme:

D 2 ∆t u n h + A h u n h - ∆t 2 12 A 2 h u n h = 0. ( 20
)
The stability condition is deduced from energy arguments and reads

∆t 2 ≤ B 00 ρ(A h ) , (21)
where B 00 = 12.

Using an Hörner algorithm, this scheme is two times more expensive than the original scheme. This over cost can be compensated by adding stabilization terms that allow to increase the time step as explained in [START_REF] Gilbert | Higher order time stepping for second order hyperbolic problems and optimal CFL conditions[END_REF] and [START_REF] Joly | Optimized higher order time discretization of second order hyperbolic problems: Construction and numerical study[END_REF].

Modified equation technique for the θ-scheme

We now use the ideas of the modified equation technique applied to the θ-scheme (instead of the leap frog scheme, i.e θ = 0). Let u h (t) be solution of (2). The truncation error of the θ-scheme (13) is:

D 2 ∆t u h (t n ) + A h {u h (t n )} θ = ∞ m=1 ∆t 2m 2 (2m + 2)! d 2 dt 2 + 2θ (2m)! A h d 2m dt 2m u h (t n ), (22)
in which we replace the time derivatives by -A h , as in (20):

D 2 ∆t u h (t n ) + A h {u h (t n )} θ + ∞ m=1 e m (θ)∆t 2m A m+1 h u h (t n ) = 0, (23)
where the coefficients e m (θ) are defined by

e m (θ) = (-1) m 2 (2m + 2)! - 2θ (2m)! . ( 24
)
To obtain a fourth order scheme, the natural idea would be to follow the modified equation procedure by keeping the first term of the series while replacing u h (t n ) with u n h . This gives

D 2 ∆t u n h + A h {u n h } θ + (θ - 1 12 )∆t 2 A 2 h u n h = 0. (25)
It is possible to study the stability of this scheme using energy techniques. As this will be a special case of the analysis done below in section 6.1, we just give here the result : the time step restriction is given by

∆t 2 ≤ B 0 (θ) ρ(A h ) , (26)
where

B 0 (θ) =          12 (1 -12θ) if θ ≤ 2 - √ 3 4 √ 3 , r(θ, 0) if θ > 2 - √ 3 4 √ 3 , (27)
where r(θ, 0) will be introduced in the next sections and satisfies r(θ, 0) 10, for θ > -(2 + √ 3)/(4 √ 3). These restrictions turn out to be quite penalizing since, unlike with the classical θ-scheme, there is no possible choice of θ that leads to an unconditionally stable scheme. Indeed, the condition (26) implies that ∆t 2 ρ(A h ) must bounded by

sup θ∈R B 0 (θ) = B 0 2 - √ 3 4 √ 3 = 6 2 - √ 3 ≃ 22.4. ( 28
)

New fourth order implicit schemes

Trying to improve the previous time step restriction, we choose to introduce a new real number ϕ ∈ R and to approximate u h (t n ) in (23) by {u n h } ϕ . This will lead us to consider a class of schemes parametrized by (θ, ϕ) ∈ R 2 , which includes scheme (25) when choosing ϕ = 0. We will call (θ, ϕ)-schemes the following schemes:

D 2 ∆t u n h + A h {u n h } θ + θ - 1 12 ∆t 2 A 2 h {u n h } ϕ = 0. ( 29
)

Inria

Next sections are dedicated to studying the numerical, stability and consistency properties of this general class of two parameters, at least fourth order, implicit schemes.

Remark 3.1. The value θ = 1/12 cancels out the new added term, retrieving (13). Therefore the stability condition has already been given by (14), which gives here:

∆t 2 ≤ 6 ρ(A h ) . ( 30
)
In the following, we will consider θ = 1/12.

Computation of the discrete solution

Scheme (29) does not seem easy to invert at first sight because it involves powers of the matrix A h which can already be difficult to invert. Yet, the computation of the solution is not a secondary question to deal with when developing a numerical scheme devoted to be actually used. We propose in this preliminary section an enhancing algorithm based on the complex factorization of a polynomial function. It will allow us at the same time to introduce notations and mathematical objects that will be widely used hereafter during convergence, energy and stability analysis.

The first step towards implementation is to rewrite the numerical scheme in a way that emphasizes its implicit nature. Assuming that ϕ = 0, and multiplying (29) by ϕ ∆t 2 , we can use the algebraic relations ( 8) and ( 9) to obtain the equivalent scheme

I h + θ∆t 2 A h + ϕ(θ - 1 12 )∆t 4 A 2 h {u n h } ϕ = u n h + [θ -ϕ]∆t 2 A h u n h =: b n h . (31)
The computational algorithm associated to the scheme (31) follows two stages:

• We retrieve {u n h } ϕ knowing u n h by inverting the matrix polynomial

P (∆t 2 A h ; θ, ϕ) = I h + θ∆t 2 A h + ϕ(θ - 1 12 )∆t 4 A 2 h . (32)
The invertibility of the matrix P (∆t 2 A h ; θ, ϕ) is equivalent to the existence of a discrete solution.

• We then use the knowledge of {u n h } ϕ to compute u n+1 h :

u n+1 h = {u n h } ϕ + (2ϕ -1)u n h ϕ -u n-1 h . ( 33
)
The first stage of this algorithm is not straightforward because the matrix to invert involves powers of the matrix A h . Indeed, if A h has a band structure (as in most classical finite elements methods), the bandwidth of A 2 h will be even larger, which penalizes the use of direct solvers. In the same way, iterative methods may suffer from the bad conditioning of P (∆t 2 A h ; θ, ϕ) since ρ(A 2 h ) = ρ(A h ) 2 (see for instance remark 7.1). To overcome such difficulties, we propose an algorithm based on the factorization of P (λ, ; θ, ϕ). Different cases arise according to the nature of its roots (r + , r -):

• The roots are real. In this case invertibility is not granted and must be verified. If P (∆t 2 A h ; θ, ϕ) is indeed invertible, {u n h } ϕ can be obtained in two steps, solving the linear systems with an intermediary unknown v n h :

v n h = (r + I h -∆t 2 A h ) -1 b n h , {u n h } ϕ = r + r -(r -I h -∆t 2 A h ) -1 v n h . (34)
• The roots are conjugate complex numbers: r + = r -. As above, we introduce v n h computed as:

v n h = (r + I h -∆t 2 A h ) -1 b n h . ( 35
)
As b n h is a real vector, we know that the solution {u n h } ϕ will also be a real vector. Therefore we can identify the imaginary parts of both sides in

1 |r + | 2 (r + I h -∆t 2 A h ){u n h } ϕ = v n h , (36)
to get, without any additional matrix inversion:

{u n h } ϕ = - |r + | 2 ℑm(r + ) ℑm(v n h ). (37)
This method only requires one matrix inversion which however happens in the complex domain. Adapted techniques can be used to avoid over cost (see [START_REF] Axelsson | Real valued iterative methods for solving complex symmetric linear systems[END_REF], [START_REF] Freund | On conjugate gradient type methods and polynomial preconditioners for a class of complex non-Hermitian matrices[END_REF]).

Remark 4.1. Using the factorization of the polynomial matrix

P (∆t 2 A h ; θ, ϕ), if A h = M -1 h K h one needs to compute the solution of linear problems of the form (r ± I h -∆t 2 M -1 h K h )x h = b h ⇐⇒ (r ± M h -∆t 2 K h )x h = M h b h . ( 38
)
This shows that the mass matrix inversion is included in the regular matrix inversion needed by the implicit scheme. This leads to nearly no over cost even if the mass matrix is non-diagonal.

Remark 4.2. To take into account the non-homogeneous version of (2), namely

d 2 dt 2 u h + A h u h = f h , u h (0) = u 0,h , du h dt (0) = u 1,h , (39)
where f h is a source term, the fully discrete (31) must be modified so that the fourth order accuracy is preserved. Using the same arguments as those used to obtained (16) from (15), one can show that, when ϕ = 0, the correct discrete scheme is

P (∆t 2 A h ; θ, ϕ){u n h } ϕ = b n h + ϕ ∆t 2 f h (t n ) + ϕ ∆t 4 12 d 2 dt 2 f h (t n ) + ϕ ∆t 4 θ - 1 12 A h f h (t n ).

Energy preservation and stability

Using relation (10), our (θ, ϕ)-schemes (29) are equivalent to:

I h + (θ - 1 4 )∆t 2 A h + (θ - 1 12 )(ϕ - 1 4 )∆t 4 A 2 h D 2 ∆t u n h + A h + (θ - 1 12 )∆t 2 A 2 h {u n h } 1/4 = 0. ( 40
)
This simplifies to the general form:

M h D 2 ∆t u n h + K h {u n h } 1/4 = 0, (41)
Inria where M h and K h are symmetric matrices defined by

K h = A h Q 2 (∆t 2 A h ; θ, ϕ), and M h = Q 1 (∆t 2 A h ; θ, ϕ) 4 , ( 42
) with      Q 1 (λ; θ, ϕ) = 4 + (4θ -1)λ + (4ϕ -1)(θ - 1 12 )λ 2 , Q 2 (λ; θ, ϕ) = 1 + (θ - 1 12 )λ. (43)
Taking the euclidian scalar product (•, •) of equation ( 41) with

(u n+1 h -u n-1 h
)/2∆t gives the energy preservation

E n+1/2 -E n-1/2 ∆t = 0, (44)
where the discrete energy E n+1/2 is defined by

E n+1/2 = 1 2 M h u n+1 h -u n h ∆t , u n+1 h -u n h ∆t + 1 2 K h u n+1 h + u n h 2 , u n+1 h + u n h 2 . ( 45
)
One can prove that the positivity of the matrices M h and K h leads to the stability of the scheme (41). Since those matrices depend on A h through the polynomials Q 1 and Q 2 , a sufficient condition can be derived to ensure this positivity:

Q 1 (λ; θ, ϕ) ≥ 0 and Q 2 (λ; θ, ϕ) ≥ 0, ∀λ ∈ [0, ∆t 2 ρ(A h )]. (46)
These conditions lead to an upper bound on ∆t 2 that depends on the values of θ = 1/12 and ϕ as stated in the following theorem:

Theorem 5.1 (CFL condition). The matrices Q 1 (∆t 2 A h ; θ, ϕ) and Q 2 (∆t 2 A h ; θ, ϕ) are positive matrices if ρ(A h ) ∆t 2 ≤ B(θ, ϕ) = min B Q1 (θ), B Q2 (θ, ϕ) (47)
where

B Q2 (θ) =      +∞ if θ > 1 12 , 12 1 -12θ otherwise, (48)
B Q1 (θ, ϕ) =            +∞ if (θ, ϕ) ∈ D US , 4 1 -4θ if (θ, ϕ) ∈ I 1/4 , r(θ, ϕ) otherwise, (49)

RR n°8090

and

I 1/4 = (-∞, 1/4) ∪ {1/4}, D US = [1/4, +∞) × [1/4, +∞) ∪ S -∪ S + , S -= ϕ ≤ 1 4 [1 + (4θ -1) 2 16(θ -1/12) ], θ < 1/12 , S + = ϕ ≥ 1 4 [1 + (4θ -1) 2 16(θ -1/12) ], θ > 1/12 , ∆(θ, ϕ) = (4θ -1) 2 -16 (4ϕ -1)(θ -1/12), r(θ, ϕ) = 1 -4θ -∆(θ, ϕ) 2(4ϕ -1)(θ -1/12
) .

(50)

Proof. The proof of this statement will be given in annex. Remark 5.1. The coefficients B 0 of (26) and B 00 of (21) can be seen as restrictions of B on specific areas of the R 2 plane:

B 0 (θ) = B(θ, 0), B 00 = B(0, 0). (51)
The following result states that respecting the CFL condition almost always leads to a wellposed discrete problem:

Theorem 5.2 (Existence of the discrete solution). If Q 1 (∆t 2 A h ; θ, ϕ) and Q 2 (∆t 2 A h ; θ, ϕ) are positive matrices, then P (∆t 2 A h ; θ, ϕ) is invertible if (θ, ϕ) = 1 12 - 1 ∆t 2 ρ(A h ) , 1 12 . ( 52
)

Inria

Proof. We write the polynomial P (λ; θ, ϕ) as a sum of two positive terms:

P (λ; θ, ϕ) = Q 1 (λ; θ, ϕ) + λQ 2 (λ; θ, ϕ) 4 . ( 53
)
For λ ∈ [0, ∆t 2 ρ(A h )], it vanishes if and only if both terms vanish at the same point. It can only

happen if λ = ∆t 2 ρ(A h ) and if (θ, ϕ) = 1/12 -(∆t 2 ρ(A h )) -1 , 1/12 .
6 Peculiar (θ, ϕ)-schemes 6.1 A class of optimal (θ, ϕ)-schemes

In the following, we are going to find the "best possible" stable (θ, ϕ)-schemes, which means to find the optimal values of θ and ϕ in the R 2 plane that minimize the consistency error of the scheme, under the constraints of scheme stability. Indeed, we know (by construction) that these schemes are fourth order accurate in time, and the consistency errors of order six and eight depend on the values of θ and ϕ. The consistency error of scheme ( 29) is obtained, first by evaluating the approximation (6) in the solution u h (t n ) of (2):

{u h (t n )} ϕ = u h (t n ) - ∞ m=1 c m (ϕ)∆t 2m A m h u h (t n ), (54)
where the coefficients c m (ϕ) are defined by

c m (ϕ) = (-1) m+1 2ϕ (2m)! . (55)
Then, the choice to approximate u h (t n ) in (23) by {u n h } ϕ leads to

D 2 ∆t u h (t n ) + A h {u h (t n )} θ + e 1 (θ)∆t 2 A 2 h {u h (t n )} ϕ + ∞ m=1 c m (ϕ)∆t 2m A m h u h (t n ) + e 2 (θ)∆t 4 A 3 h u h (t n ) + e 3 (θ)∆t 6 A 4 h u h (t n ) = O(∆t 8 ), (56)
which gives

D 2 ∆t u h (t n )+A h {u h (t n )} θ +e 1 (θ)∆t 2 A 2 h {u h (t n )} ϕ = -ε 3 (θ, ϕ)∆t 4 A 3 h u h (t n )-ε 4 (θ, ϕ)∆t 6 A 4 h u h (t n )+O(∆t 8 ), (57)
where the first terms of the consistency error are

         ε 3 (θ, ϕ) = e 2 (θ) + e 1 (θ) c 1 (ϕ) = 1 360 - θ 12 - ϕ 12 + θ ϕ, (58a)
ε 4 (θ, ϕ) = e 3 (θ) + e 1 (θ) c 2 (ϕ) = -1 20160 + θ 360 + ϕ 144 - θ ϕ 12 . ( 58b
)
Regarding the stability conditions, we will see that they provide nonlinear constraints on (θ, ϕ) which depend on ∆t 2 ρ(A h ). This is why we tackle this issue using a non standard point of view : we assume ∆t 2 ρ(A h ) to be known and we solve the corresponding optimization problem:

min (θ,ϕ) ∈ R 2 |ε 3 (θ, ϕ)|, Q 1 (λ; θ, ϕ) ≥ 0, Q 2 (λ; θ, ϕ) ≥ 0, ∀λ ∈ [0, ∆t 2 ρ(A h )]. ( 59
)
If it is possible to find values of (θ, ϕ) that make ε 3 vanish in the stability region, then the optimal choices within these values are the one that minimize the absolute value of ε 4 .

In section 6.1.1, we will try to mimic the super-convergence phenomenon that we find in the classical θ-scheme when θ = 1/12. Indeed, this second order scheme appears to be fourth order accurate for the peculiar choice of θ = 1/12. In our case, we will see that it will be possible to obtain stable sixth order schemes by restraining the choice of the couple (θ, ϕ) to a curve in the R 2 plane corresponding to the zeros of (58a), and even stable eighth order schemes by choosing a special couple on this curve that correspond to a zero of (58b). The major drawback of this powerful result is that it can only happen for small values of ∆t 2 ρ(A h ).

When this approach is not possible (for ∆t 2 ρ(A h ) greater than a certain value), section 6.1.2 investigates which couple (θ, ϕ) in the R 2 plane leads to the stable numerical scheme that minimizes the absolute value of the consistency error (58a).

Finally in section 6.1.3, we construct "optimal" unconditionally stable schemes by assuming in the optimization process that ∆t 2 ρ(A h ) = +∞. This will lead to (θ, ϕ)-schemes, where θ and ϕ depend on a small parameter δ, and that minimize the absolute value of (58a) when δ tends to 0.

Remark 6.1. It is possible to obtain higher order θ-schemes parametrized by more real numbers.

For instance, if we extend the previous approach to the next order of approximation, we can introduce (θ, ϕ, ψ, ψ) ∈ R 4 and the obtained schemes read:

D 2 ∆t u n h + A h {u n h } θ + e 1 (θ)∆t 2 A 2 h {u n h } ϕ + ∆t 4 A 3 h e 2 (θ){u n h } ψ + e 1 (θ)c 1 (ϕ){u n h } ψ = 0. ( 60
)
The consistency error of these schemes is given by

-∆t 6 A 4 h e 3 (θ) + e 1 (θ)c 1 (ϕ)c 1 ( ψ) + e 1 (θ)c 2 (ϕ) + e 2 (θ)c 1 (ψ) u h (t n ) + O(∆t 8 ), (61)
showing that these schemes are at least sixth order accurate.

6.1.1 Sixth and eighth order stable schemes.

We look for stable (θ, ϕ)-schemes such that ε 3 vanishes and such that |ε 4 | is minimized. This restricts the choice of θ and ϕ to a curve in R 2 described by

ϕ ⋆ = (12 θ ⋆ -1) -1 (θ ⋆ - 1 30 ). (62)
All the values on this curve lead to a sixth order scheme. The stability conditions require that for all λ ∈ [0, ∆t 2 ρ(A h )] we must have

       Q 1 (λ; θ ⋆ , ϕ ⋆ ) = (4θ ⋆ -1)λ + (- 2 3 θ ⋆ + 13 180 )λ 2 + 4 ≥ 0, Q 2 (λ; θ ⋆ , ϕ ⋆ ) = 1 + (θ ⋆ - 1 12 )λ ≥ 0, (63)
where ϕ has been eliminated using (62). The second inequality is constraining only if θ ⋆ < 1/12, in which case

∆t 2 ρ(A h ) ≤ 12 1 -12θ ⋆ := p(θ ⋆ ). (64)

Inria

The first inequality is respected if and only if the polynomial Q 1 (λ; θ ⋆ , ϕ ⋆ ) has no root of multiplicity one in the interval [0, ∆t 2 ρ(A h )]. Let ∆(θ ⋆ ) be the discriminant of Q 1 (λ; θ ⋆ , ϕ ⋆ ):

∆(θ ⋆ ) = 16 (θ ⋆ ) 2 + 8 3 θ ⋆ - 7 45 . ( 65
)
Different cases arise according to the convexity of Q 1 :

• Concave when θ ⋆ > 13 120 . In this case, there are two simple roots of opposite signs. The positive root must be greater than ∆t 2 ρ(A h ), that is:

∆t 2 ρ(A h ) ≤ r(θ ⋆ ) := 1 -4θ ⋆ -∆(θ ⋆ ) 2 (13/180 -2θ ⋆ /3) . (66)
• Linear when θ ⋆ 0 = 13 120 . In this case, Q 1 is a decreasing linear function, which is positive on [0, 120 17 ] to which ∆t 2 ρ(A h ) must belong. This condition extends by continuity relation (66) to its singular point.

• Convex when θ ⋆ < 13 120 . In this case, we have Q ′ 1 (0; θ ⋆ , ϕ ⋆ ) < 0 and two different cases according to the sign of ∆(θ ⋆ ) must be considered:

• ∆(θ ⋆ ) > 0 ⇔ θ ⋆ / ∈ [-1 12 - √ 15 30 , -1 12 + √ 15
30 ] : Two simple roots. Therefore the lower root must be greater than ∆t 2 ρ(A h ). This condition turns out to be the same as (66). Figure 2: Bounds on ∆t 2 ρ(A h ) for any θ ⋆ : the scheme is stable if the couple (∆t 2 ρ(A h ), θ ⋆ ) lies in both hatched areas. In continuous line, we represented p(θ ⋆ ) while in dotted line we represented r(θ ⋆ ) when ∆(θ ⋆ ) > 0.

• ∆(θ ⋆ ) ≤ 0 ⇔ θ ⋆ ∈ [-1 12 - √ 15
To sum up, the upper bound on ∆t 2 ρ(A h ) is :

                       min(p(θ ⋆ ), r(θ ⋆ )) if θ ⋆ < - 1 12 - √ 15 30 , p(θ ⋆ ) if - 1 12 - √ 15 30 ≤ θ ⋆ ≤ - 1 12 + √ 15 30 , min(p(θ ⋆ ), r(θ ⋆ )) if - 1 12 + √ 15 30 < θ ⋆ < 1 12 , r(θ ⋆ ) if 1 12 < θ ⋆ . (67)
In this variety of sixth order schemes, the choice of θ ⋆ is guided by the minimization of the absolute value of ε 4 , which reads

ε 4 (θ ⋆ , ϕ ⋆ ) = -θ ⋆ 240 + 11 60480 . ( 68
)
It is a linear decreasing function of θ ⋆ , vanishing for the value

θ ⋆⋆ = 11 252 ⇒ ϕ ⋆⋆ = - 13 600 (69)
• For ∆t 2 ρ ∈ [0, 126 5 ], the couple (∆t 2 ρ(A h ), θ ⋆⋆ ) provides an eighth order accurate stable scheme.

• For ∆t 2 ρ ∈ [ 126 5 , 30

5-√

15 ], the couple (∆t 2 ρ(A h ), θ ⋆⋆ ) violates the constraints (67). We can prove that minimizing |ε 4 (θ ⋆ )| while satisfying (67) leads to the choice :

θ ⋆ = 1 12 - 1 ∆t 2 ρ(A h ) . (70)
• For ∆t 2 ρ > 30 5-√ 15 , it is not possible to construct sixth order stable schemes anymore. Indeed, all the upper bounds of (67) are lower than 30 5-√ 15 .

Fourth order optimal stable scheme

We consider now the values of ∆t 2 ρ(A h ) which could not satisfy the stability conditions of previous section, more precisely we assume that ∆t 2 ρ(A h ) > 30/(5 -√ 15). As seen before it is not possible to make ε 3 vanish in this case, thus we will construct a stable scheme for a given ∆t 2 ρ(A h ) that minimizes its absolute value by solving (59).

The second inequality of the constraints in (59) leads to

θ ≥ 1 12 - 1 ∆t 2 ρ(A h ) . ( 71
)
For the first inequality, different situations arise according to the convexity and the slope at the origin of Q 1 (λ; θ, ϕ). It divides the R 2 plane into six areas delimited by the values 1/12 and 1/4 for θ and 1/4 for ϕ.

Let us first consider the area 1/12 -1/(∆t 2 ρ(A h )) ≤ θ < 1/12, ϕ < 1/4. Q 1 is then a convex parabola with a negative slope at the origin.

Inria

• The parabola crosses the y-axis when ∆(Q 1 ) > 0. In this case, the interesting interval [0, ∆t 2 ρ(A h )] must be before the first root, or in other terms, we must have

Q 1 (∆t 2 ρ(A h )) ≥ 0 and Q ′ 1 (∆t 2 ρ(A h )) ≤ 0.
But the calculation shows that these two conditions are incompatible because we made the assumption on ∆t 2 ρ(A h ) > 8.

• The parabola never crosses the y-axis when ∆(Q 1 ) ≤ 0 which leads to the condition:

ϕ ≤ 1 4 1 + (4θ -1) 2 16 (θ -1/12) . ( 72
)
In this part of the plane, the sixth order consistency error is a decreasing function of ϕ, for any θ:

ε 3 (θ, ϕ) = 1 360 - θ 12 + ϕ (θ - 1 12 ) ( 73
)
and stays positive when (72) is respected. Hence, if θ is fixed, the value of ϕ that minimizes the absolute value of ε 3 is the greatest value possible. As a consequence, we want to minimize |ε 3 | on the curve

θ, ϕ = 1 4 1 + (4θ -1) 2 16 (θ -1/12) , ( 74
)
which leads to the 1D minimization of:

ε 3 θ, 1 4 1 + (4θ -1) 2 16(θ -1/12) = θ 2 4 + θ 24 - 7 2880 (75)
This objective function is positive and increasing on the interval [1/12 -1/∆t 2 ρ(A h ), 1/12] as soon as ∆t 2 ρ(A h ) ≥ 4. Therefore, the optimal values of θ and ϕ in this area of R 2 are given by:

θ ♯ = 1 12 - 1 ∆t 2 ρ(A h ) , ϕ ♯ = 1 4 - ∆t 2 ρ(A h ) 64 2 3 + 4 ∆t 2 ρ(A h ) 2 . ( 76
)
It turns out that in the other areas, either it is not possible to construct stable schemes, or stable schemes give a greater consistency error. We provide in annex a proof of this statement. Notice that if the product ∆t 2 ρ(A h ) gets very large, the optimal values (θ ♯ , ϕ ♯ ) tend to (1/12, -∞).

6.1.3 Fourth order unconditionally stable scheme.

When the spectral radius of the operator is not known, we assume it is infinite, which means that the positivity of Q 1 and Q 2 must be ensured on the whole R + interval. Unfortunately it is not possible to pass to the limit in the formulae (θ ♯ , ϕ ♯ ). We first notice that it is necessary to have

θ > 1 12 , ϕ > 1 4 , ( 77
)
otherwise Q 2 is a decreasing affine function and Q 1 a concave parabola. Again, we have to distinguish several cases depending on the slope at the origin of Q 1 .

• Either θ ≥ 1/4 and stability is acquired without any other condition.

• Either 1/12 < θ < 1/4 and the parabola Q 1 must not cross the y-axis : its discriminant must be negative, leading to

ϕ ≥ 1 4 1 + (4θ -1) 2 16 (θ -1/12) . ( 78
)
It both cases, ε 3 is a positive decreasing function of ϕ, therefore the lowest possible value of ϕ must be chosen. Again, we tackle a 1D optimization problem. Let us introduce δ > 0 such that

θ δ = 1 12 + δ ⇒ ϕ δ = (-2/3 + 4δ) 2 64 δ + 1 4 . ( 79
)
The consistency error parametrized by δ is given by

ε 3 (θ δ , ϕ δ ) = δ 2 4 + δ 12 + 1 360 . ( 80
)
Therefore δ must be chosen as low as possible, even if the lower bound 0 is not achievable. To sum up, we have constructed a family of unconditionally stable (θ, ϕ)-schemes parametrized by δ > 0 that minimize the consistency error "at the limit".

Conclusions

The results of previous sections are summarized in table 6.1.4. With these choices of (θ, ϕ)scheme, the polynomial P (∆t 2 ρ(A h ); θ, ϕ) defined by (32) has complex roots as soon as

∆t 2 ρ(A h ) ≥ 60 7 + 24 7 √ 15 ≃ 21.8502, (81)
otherwise the roots are negative. In both cases, the invertibility of the matrix P (∆t 2 A h , θ, ϕ) is guaranteed.

Table 1: Optimal values of θ and ϕ for a given positive ρ := ∆t 2 ρ(A h ) and for a positive small parameter δ.

"Multiple roots" fourth order schemes

The schemes developed in section 6.1 lead to a complex factorization of the polynomial (32) as soon as ∆t 2 ρ(A h ) is big enough. In this section we build unconditionally stable schemes for which the polynomial P (λ; θ, ϕ) has a double negative root r:

P (∆t 2 A h ; θ, ϕ) = r -2 (rI h -∆t 2 A h ) 2 . ( 82
)
Inria Therefore, the inversion of P (∆t 2 A h ; θ, ϕ) can be done by inverting the same matrix (rI h -∆t 2 A h ) twice, which reduces cost when direct solvers are used. This happens when the discriminant of P (λ; θ, ϕ) vanishes, giving the following relation between θ and ϕ:

ϕ = 3θ 2 12θ -1 . (83)
For this choice of ϕ, the stability conditions read

   1 + (θ -1/12)λ ≥ 0, 4 + (4θ -1)λ + (θ 2 -θ + 1 12 )λ 2 ≥ 0, ∀λ ≥ 0. (84)
The first condition implies that θ > 1/12. The second condition cannot be respected unless the leading coefficient is positive, which leads to

θ / ∈ [1/2 -1/ √ 6, 1/2 + 1/ √ 6].
Since 1/12 lies in this interval, θ must be greater than 1/2 + 1/ √ 6, implying that the negative multiple root of P (λ; θ, ϕ) is given by r = -2/θ. To sum up, the (θ, ϕ)-schemes satisfying (83) and

θ ≥ 1 2 + 1 √ 6 ≃ 0.91 (85)
are unconditionally stable and we get the factorization

P (∆t 2 A h ; θ, ϕ) = (I h + θ ∆t 2 2 A h ) 2 . (86)
For these choices of (θ, ϕ) the absolute value of the consistency error |ε 3 | is a positive and increasing function of θ, so the value θ = 1/2 + 1/ √ 6 (which corresponds to ϕ = 1/4) is optimal.

7 Numerical results

1D convergence test cases

In this section we present an academical test case for which we know the exact solution, allowing us to compute exactly numerical errors. We consider the 1D scalar wave equation with velocity 1 in a domain of length 1 with periodic boundary conditions:

∂ 2 ∂t 2 u(x, t) - ∂ 2 ∂x 2 u(x, t) = 0, x ∈ [0, 1], u(0, t) = u(1, t), ∂ ∂x u(0, t) = ∂ ∂x u(1, t). (87)
The initial condition is a gaussian pulse centered in the domain. The initial velocity is such that the analytical solution is this gaussian pulse traveling from left to right. More precisely, the initial conditions are

u(x, 0) = u 0 (x) = e (ε-1/2) -2 ln(ε) (x-1/2) 2 if |x -1/2| ≤ 1/2 -ε ε if |x -1/2| > 1/2 -ε , ∂ ∂t u(x, 0) = - ∂ ∂x u 0 (x),
where ε is a small parameter (in practice we choose ε = 10 -20 ). The analytic solution is given by u

(x, t) = u 0 (x -t + n) where n ∈ N such that x -t + n ∈ [0, 1].
In order to solve numerically this wave equation, the natural idea is to use the second order centered finite difference scheme. It is even known to provide the exact solution on regular meshes if the time step is chosen equal to its maximal value allowed by the CFL condition. In this context (regular mesh) it is obvious that our implicit schemes would not compare well. They were originally designed to overcome specific difficulties as distorted meshes or rapidly varying coefficients. In order to pedagogically represent these kinds of difficulties, we propose to introduce a very small element in the mesh. It can seem artificial in this simple 1D example, but is numerically representative of what can happen in 2D or 3D when realistic domains are involved.

Standard explicit and θ-schemes are compared with some of the new (θ, ϕ)-schemes introduced above, and with the exact solution. The error e is computed as the sup. over time of the relative discrete L 2 error. The estimated cost corresponds to the total number of matrix-vector products with A h required for the simulation. Non-preconditioned iterative methods (Conjugate Gradient and Minimal Residual, see [START_REF] Freund | On conjugate gradient type methods and polynomial preconditioners for a class of complex non-Hermitian matrices[END_REF]) have been used to solve the implicit schemes in factorized form as presented in 4.

The mesh is composed of 12 identical elements surrounding a small element of size τ = 10 -4 , which leads to a spectral radius ρ(A h ) = 5.62 × 10 9 when sixth order spectral finite elements with mass lumping are used to tackle spatial discretization (see [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF]). For the explicit scheme, the time step is restricted to 2.7 × 10 -5 by the CFL condition, whereas for the implicit schemes, we impose ∆t = 0.016 or ∆t = 0.004. The simulation runs until T = 10. Table 2 summarizes the performances of the chosen schemes and figure 3 shows the obtained snapshots at final time. Table 2: Comparison between several schemes: leap-frog explicit (θ = 0), θ-scheme with θ = 1/4, naive (θ, ϕ)-scheme with (θ, ϕ) = (1/4, 1/4) and optimal (θ, ϕ)-scheme adapted to the product ∆t 2 ρ(A h ) (given by equation ( 76)). As expected, the small element clearly penalizes the explicit scheme: a small step step must be chosen, which increases the cost of the method but also increases its accuracy. This is why Inria the final snapshot is very close to the analytical solution, and the error is very low, for a very expensive cost. For the three following implicit schemes, we choose a time step of 1.6 × 10 -2 . The θ-scheme with θ = 1/4 is very cheap but behaves very poorly, as illustrated by the snapshot and the relative error of about 0.2. The naive (θ, ϕ)-scheme with (θ, ϕ) = (1/4, 1/4) gives a much better error and a nicer snapshot, even if numerical dispersion is still visible. The optimal (θ, ϕ)-scheme clearly overcomes the naive (θ, ϕ)-scheme, by being more accurate and cheaper. The first observation was to be expected since the criterion of optimality was indeed accuracy, but the second one had not been foreseen and is linked to the correlation between the condition number of the matrix P (∆t 2 A h ; θ, ϕ) and the consistency error ε 3 . Finally, by choosing an appropriate time step of 4.0 × 10 -3 , we recover with the optimal (θ, ϕ)-scheme the accuracy of the explicit scheme, while staying about five times cheaper.

Remark 7.1. To highlight the interest of factorizing the polynom P in this specific 1D case, let us focus on the (θ, ϕ)-scheme with (θ, ϕ) = (1/4, 1/4). On the one hand, the natural approach would be to invert the matrix

P (∆t 2 A h ; 1/4, 1/4) = I h + ∆t 2 4 A h + ∆t 4 24 A 2 h .
Its condition number is given by P (∆t 2 ρ(A h ); 1/4, 1/4) = 8.62 × 10 10 . On the other hand, we can use the factorized form introduced in section 4. The complex roots of the polynomial P read r ± = -3 ± i √ 15, hence the matrix to invert is (-3 + i √ 15)I h -∆t 2 A h and its condition number is equal to |r + -∆t 2 ρ(A h )|/|r + | = 2.94 × 10 5 . This last matrix is far easier to invert both in terms of computational cost (iterative methods) and numerical precision.

A convergence study has been lead on a regular mesh with sixth order finite elements (or seventh order to obtain an eighth order convergence curve) in space for different (θ, ϕ)-schemes. The error is computed as the sup norm in time (here T = 8) of the discrete L 2 norm (induced by the mass matrix) of difference with the analytical solution. The results are given in figure 4.

Figure 4: Log of the error w.r.t the log of the mesh size for different (θ, ϕ)-schemes and discretization parameters. + : optimal "multiple roots" scheme with ∆t 2 ρ = 120, • : optimal (θ, ϕ)-scheme with ∆t 2 ρ = 120, : optimal (θ, ϕ)-scheme with ∆t 2 ρ = 60, × : optimal (θ, ϕ)-scheme with ∆t 2 ρ = 26, • : optimal (θ, ϕ)-scheme with ∆t 2 ρ = 22.

We obtain the expected rates of convergence for the different (θ, ϕ)-schemes. All the simulations and matrix inversions have been done with an iterative method without preconditioning. This is clearly non-optimal in term of computational cost but also in term of accuracy since it induced a numerical locking that prevented us to reach better accuracy with the eighth order scheme.

Propagation of acoustic waves in a strongly heterogeneous 2D domain

To conclude the numerical results, we present a more complex configuration. We are interested in the propagation of acoustic pressure p(x, t) in a square domain Ω = [-1, 1] × [-1, 1] in which we assume that the propagation is governed by the non-homogeneous equations

∂ 2 ∂t 2 p(x, t)-div c(x) ∇p(x, t) = f (x, t), x ∈ Ω, ∇p•n = 0, x ∈ ∂Ω, ∂ ∂t p(x, 0) = p(x, 0) = 0,
where n is the outward normal, c(x) the material properties of the domain of propagation and f (x, t) := g(x) h(t) a source term. The parameters c(x) and h(t) are represented in figure 5 whereas g(x) is a 2D gaussian pulse centered around x = 0. The explicit expressions of these coefficients are c(x) = 1+10 4 e -800 (|x| 2 -0.35) 2 , h(t) = 2π 2 (16 (t-0.1)-1) 2 -1 e -π 2 (16 (t-0.1)-1) 2 , g(x) = e -800 |x| 2 . Explicit computation of the discrete solution of this problem may not be adapted if one wants to minimize the computational cost while reaching a target, but not necessarily very small, accuracy. Indeed, the upper bound for the time step in explicit methods is proportional to 1/ max x c(x) = 10 -2 . This is a typical situation where implicit and unconditionally stables schemes are often used, since in this case the time step is no longer related to c(x) and can be tuned to reach the desired accuracy and obtain a reasonable computational cost.

We choose to compare three different types of discretization in time to a precomputed refined reference solution:

• Schemes A & B Two standard leap frog schemes, of second and fourth order, which are explicit. Hence they must be used with a very small time step. They should yield a very good accuracy yet being more costly than the implicit schemes. Note that the fourth order Inria leap frog scheme corresponds to the choice (θ, ϕ) = (0, 0) also called modified equation scheme.

• Scheme C The second order θ-scheme with θ = 1/4. It is the most classical scheme that can be used to tackle the implicit discretization of our problem.

• Scheme D Our fourth order discretization (θ, ϕ)-scheme with (θ, ϕ) = (1/2 + 1/ √ 6, 1/4) as computed in subsection 6.2. This is the most easily implantable alternative that we presented and it should outperform the second-order θ-scheme.

The same space discretization is chosen for all the numerical computations since we are only interested in comparing the time discretization of a semi-discrete problem similar to equation (39). The spatial mesh is a fine enough uniform grid of size h = 0.025 (1600 squared elements) and we use fifth order spectral finite elements with mass lumping (see [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF]) (160801 degrees of freedom and ρ(A h ) ≃ 1.2 × 10 10 ). Note also that the source term is handled as indicated in remark 4.2 to avoid any accuracy loss.

The explicit schemes A and B must respect their stability condition, leading to a maximal time step allowed of respectively ∆t ≈ 1.8 × 10 -5 and ∆t ≈ √ 3 × 1.8 × 10 -5 . The last two schemes C and D are unconditionally stable. The time step is chosen to be ∆t = 1.8 × 10 -3 , which is close to the maximal time step allowed by the explicit scheme if the coefficient c(x) was constant (equal to 1) over all of the domain. These schemes require one (respectively two) inversions of a real matrix at each time step. The choice of the peculiar fourth order scheme D described in subsection 6.2 has been motivated by the fact that the same (symmetric positive definite) matrix must be inverted twice at each time step, which implies that in practice only one LU factorization is performed before the iterations start. Some snapshots of the solution are presented in figure 6 and the results of the numerical experiments are summed-up in table 3. The reference solution is computed with the fourth order classical θ-scheme B, with a time step of ∆t = 3.11 × 10 -6 , which is ten times less than the maximum value allowed. Since the time steps taken for the explicit schemes A and B are very low, we see that as expected, they give a very accurate solution for a higher computational cost. More importantly we see that the implicit θ-schemes C give a solution accurate from around 9% to 2% for a cost around 7 to 4 time less, however for nearly the same computational cost (θ, ϕ)-schemes D give solutions accurate from around 4% to 0.3%, outperforming the classical θ-scheme.

Conclusions and prospects

We have introduced a new family of fourth order accurate in time, implicit, energy preserving schemes, denoted as the (θ, ϕ)-schemes, that generalize the famous second order accurate θschemes. Their stability properties have been investigated via an energy analysis, leading to a general CFL condition. We have also provided in this family schemes that minimize the consistency error. We have found schemes with super-convergence properties that can be used if ∆t 2 ρ(A h ) is small enough, or optimal schemes if ∆t 2 ρ(A h ) is large. We have also provided unconditionally stable schemes and finally schemes that require the same numerical treatment as the usual θ-scheme but are fourth order accurate. Numerical results show the interest of such schemes as well as their convergence behavior.

However, we have sometimes noted the need to use preconditioning strategies when a direct factorization of the finite elements matrices was not possible. Several authors have already suggested solutions to deal with the inversion of the complex linear system occurring with the optimal (θ, ϕ)-scheme ([AK00], [START_REF] Freund | On conjugate gradient type methods and polynomial preconditioners for a class of complex non-Hermitian matrices[END_REF], [START_REF] Benzi | Block preconditioning of real-valued iterative algorithms for complex linear systems[END_REF]). These approaches still need to be tested and adapted to the (θ, ϕ)-scheme. A natural extension of this work would be to consider and analyze the sixth order accurate schemes (see remark 6.1), which offer even more parameters and are therefore expected to give a very low consistency error after the optimization process. Finally, significant extension of this work will be to take into account dissipative terms and non trivial boundary conditions in the equations, which will be the subject of forthcoming work.

Inria

Figure 1 :

 1 Figure 1: Graphical representation of the CFL condition B(θ, ϕ). The dotted line represents the demarcation of S + . The white areas stand for an infinite upper bound (unconditionally stable schemes) whereas the colored scale stands for values of B (the darker being the lower).

  No simple root. Therefore the first condition is automatically fulfilled.

  Scheme∆t Error Cost Explicit 2.7 × 10 -5 2.4 × 10 -5 374 957 θ = 1/4 1.6 × 10 -2 2.3 × 10 -1 37 158 (θ, ϕ) = (1/4, 1/4) 1.6 × 10 -2 2.8 × 10 -2 67 526 (θ, ϕ) = (θ ♯ , ϕ ♯ ) 1.6 × 10 -2 4.1 × 10 -3 47 004 (θ, ϕ) = (θ ♯ , ϕ ♯ ) 4.0 × 10 -3 3.1 × 10 -5 73 214

Figure 3 :

 3 Figure 3: Snapshots of the numerical solutions at final time.Figure (a) corresponds to the explicit scheme, (b) to the θ-scheme, (c) to the naive (θ, ϕ)-scheme and (d) to the optimal (θ, ϕ)-scheme. All implicit schemes use the time step ∆t = 1.6 × 10 -2 .

  Figure 3: Snapshots of the numerical solutions at final time.Figure (a) corresponds to the explicit scheme, (b) to the θ-scheme, (c) to the naive (θ, ϕ)-scheme and (d) to the optimal (θ, ϕ)-scheme. All implicit schemes use the time step ∆t = 1.6 × 10 -2 .

Figure 5 :

 5 Figure 5: Representation of the coefficients c(x) and h(t).

Figure 6 :Table 3 :

 63 Figure 6: From left to right, top to bottom: snapshots of the pressure field computed with the explicit scheme A, for T = 0.105, T = 0.3, T = 0.5 and T = 0.7, T = 0.9, T = 1.1. The dashed circle is parametrized by |x| 2 = 0.35 and represents the region where c(x) reaches his maximum value.

RR n°8090

Avenue de la Vieille Tour, 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

A Proof of theorem (5.1)

The positivity of the energy is granted as soon as Q 1 (λ; θ, ϕ) and Q 2 (λ; θ, ϕ) are positive for any λ ∈ [0, ∆t 2 ρ(A h )]:

The positivity of each polynomial leads to upper bounds on ∆t 2 ρ(A h ), therefore the positivity of the energy is fulfilled if ∆t 2 ρ(A h ) is lower than the minimum of both bounds.

Positivity of Q 2 This polynomial is an affine function, being 1 at the origin. If its leading coefficient is positive (if θ > 1/12), then (88b) holds. In the opposite case, (88b) holds if and only if Q 2 (∆t 2 ρ(A h ); θ, ϕ) ≥ 0. These results lead to the first upper bound (48).

Positivity of Q 1 This polynomial is an affine function if ϕ = 1/4. In this case, either θ > 1/12 and (88a) is always true, either θ < 1/12 and (88a) holds if and only if Q 1 (∆t 2 ρ(A h ); θ, ϕ) ≥ 0, which leads to the upper bound 4/(1 -4θ). Assume now that ϕ = 1/4. The polynomial Q 1 being 4 at the origin, several possible situations arise according to the sign of the leading coefficient and the sign of the slope at the origin:

The three concave situations lead to the following upper bound: ∆t 2 ρ(A h ) must be lower than the positive root of Q 1 . The formula for this positive root r(θ, ϕ) is given by (50). In the case convex with positive slope, Q 1 is then positive at the origin and increasing, thus the condition (88a) is automatically fulfilled, therefore we set the upper bound to +∞. Finally, the two convex with negative slope situations are a little more complicated. Two cases arise according to the sign of the discriminant ∆(θ, ϕ) defined in (50). Either there are no distinct roots when ∆(θ, ϕ) ≤ 0 (if (θ, ϕ) ∈ S + or (θ, ϕ) ∈ S -), thus (88a) is fulfilled and we set the upper bound to +∞. Or there are two positive roots when ∆(θ, ϕ) > 0, therefore the condition (88a) is fulfilled if ∆t 2 ρ(A h ) is lower than the first root of Q 1 , which is still r(θ, ϕ) since the sign of the leading coefficient has changed.

B Optimal order 4 stable scheme

In this appendix, we complete the proof given in 6.1.2 which states that the optimal value of |ε 3 (θ, ϕ|)

RR n°8090

It is shown in 6.1.2 that in this quadrant, this choice gives the optimal value of |ε 3 |:

which is lower than 1/360 because ∆t 2 ρ(A h ) > 30/(5 -√ 15). We will see that in the other regions of the (θ, ϕ)-plane, either the schemes cannot be stable, either the stable schemes give a greater consistency error.

The quadrant {θ < 1/12, ϕ ≥ 1/4} The positivity of Q 2 on the interval [0, ∆t 2 ρ(A h )] imposes that θ > 1/12 -1/∆t 2 ρ(A h ). Moreover, the positivity of the concave polynomial Q 1 on the same interval is respected if Q 1 (∆t 2 ρ(A h ), θ, ϕ) ≥ 0. These two conditions are incompatible because ∆t 2 ρ(A h ) > 30/(5 -√ 15). Therefore, there are no stable schemes in this quadrant.

The half space {θ > 1/12} The positivity of Q 2 is granted in this area. As for Q 1 , we divide the space into three zones:

• ϕ < 1/4: in this area, Q 1 (λ; θ, ϕ) is a concave parabola whose positivity on the interval [0, ∆t 2 ρ(A h )] is acquired if and only if Q 1 (∆t 2 ρ(A h ); θ, ϕ) ≥ 0, which gives the following inequality:

which can be respected only if θ > 1/4 -1/∆t 2 ρ(A h ).

• ϕ ≥ 1/4 and θ ≥ 1/4: in this area, Q 1 (λ; θ, ϕ) is a convex parabola with a positive slope at the origin, being 4 at the origin. It is then positive for any λ ≥ 0. All schemes are stable in this region.

• ϕ ≥ 1/4 and 1/12 < θ < 1/4: in this area, Q 1 (λ; θ, ϕ) is a convex parabola with a negative slope at the origin. Its positivity on [0, ∆t 2 ρ(A h )] can be acquired either if there are no distinct roots (∆(θ, ϕ) ≤ 0) or if the first root r(θ, ϕ) is greater than ∆t 2 ρ(A h ). The second condition can be written in a way that avoids the inversion of the relation r(θ, ϕ):

Inria

It is easy to show that

Therefore, the lower bound on ϕ will be ϕ

To sum up, the stable schemes of the half plane θ > 1/12 are obtained in the region

Since the zero level set of ε 3 (θ, ϕ) is outside this stability zone, the minimum value is achieved on the boundary of the zone, which leads to the following optimization problem:

.

(98) This function is continuous and increasing and can be extended by continuity up to θ = 1/12 : the limit value is 1/360, which is always greater than the optimal value found in the other quadrant. This concludes our proof.


Nicolas Ballas

Benjamin Labbé

Aymen Shabou

Hervé Le Borgne

Philippe-Henri Gosselin

Miriam Redi

Bernard Merialdo

Hervé Jégou

Jonathan Delhumeau

Rémi Vieux

Bernard Mérialdo

Boris Mansencal

Jenny Benois-Pineau

Stéphane Ayache

Abdelkader Haadi

Bahjat Safadi

Franck Thollard

Nadia Derbas

Georges Quénot

Hervé Bredin

Matthieu Cord

Boyang Gao

Chao Zhu

Yuxing Tang

Emmanuel Dellandreav

Charles-Edmond Bichot

Liming Chen

Alexandre Benoît

Patrick Lambert

Tiberius Strat

Joseph Razik

Sébastion Paris

Hervé Glotin

Tran Ngoc Trung

Dijana Petrovska

Géerard Chollet

Andrei Stoian

Michel Crucianu

IRIM at TRECVID 2012: Semantic Indexing and Instance Search

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

Abstract

The IRIM group is a consortium of French teams working on Multimedia Indexing and Retrieval. This paper describes its participation to the TRECVID 2012 semantic indexing and instance search tasks. For the semantic indexing task, our approach uses a six-stages processing pipelines for computing scores for the likelihood of a video shot to contain a target concept. These scores are then used for producing a ranked list of images or shots that are the most likely to contain the target concept. The pipeline is composed of the following steps: descriptor extraction, descriptor optimization, classification, fusion of descriptor variants, higher-level fusion, and re-ranking. We evaluated a number of different descriptors and tried different fusion strategies. The best IRIM run has a Mean Inferred Average Precision of 0.2378, which ranked us 4th out of 16 participants.

For the instance search task, our approach uses two steps. First individual methods of participants are used to compute similrity between an example image of instance and keyframes of a video clip. Then a two-step fusion method is used to combine these individual results and obtain a score for the likelihood of an instance to appear in a video clip. These scores are used to obtain a ranked list of clips the most likely to contain the queried instance. The best IRIM run has a MAP of 0.1192, which ranked us 29th on 79 fully automatic runs.

1 Semantic Indexing

Introduction

The TRECVID 2012 semantic indexing task is described in the TRECVID 2012 overview paper [START_REF] Smeaton | Evaluation campaigns and TRECVid[END_REF][START_REF] Over | TRECVID 2012 -An Overview of the Goals, Tasks, Data, Eval[END_REF]. Automatic assignment of semantic tags representing high-level features or concepts to video segments can be fundamental technology for filtering, categorization, browsing, search, and other video exploitation. New technical issues to be addressed include methods needed/possible as collection size and diversity increase, when the number of features increases, and when features are related by an ontology. The task is defined as follows: "Given the test collection, master shot reference, and concept/feature definitions, return for each feature a list of at most 2000 shot IDs from the test collection ranked according to the possibility of detecting the feature." 346 concepts have been selected for the TRECVID 2012 semantic indexing task. Annotations on the development part of the collections were provided in the context of a collaborative annotation effort [START_REF] Safadi | Quaero at TRECVID 2011: Semantic Indexing and Collaborative Annotation[END_REF].

Fifteen French groups (CEA-LIST, CNAM, ETIS, EURECOM, INRIA, LABRI, LIF, LIG, LIMSI, LIP6, LIRIS, LISTIC, LSIS, Mines Télécom and Télécom ParisTech) collaborated to participate to the TRECVID 2012 semantic indexing task.

The IRIM approach uses a six-stages processing pipeline that compute scores reflecting the likelihood of a video shot to contain a target concept. These scores are then used for producing a ranked list of images or shots that are the most likely to contain the target concept. The pipeline is composed of the following steps:

1. Descriptor extraction. A variety of audio, image and motion descriptors have been produced by the participants (section 1.2).

2. Descriptor optimization. A post-processing of the descriptors allows to simultaneaously improve their performance and to reduce their size (section 1.3).

3. Classification. Two types of classifiers are used as well as their fusion (section 1.4).

4. Fusion of descriptor variants. We fuse here variations of the same descriptor, e.g. bag of word histograms with different sizes or associated to different image decompositions (section 1.6).

5. Higher-level fusion. We fuse here descriptors of different types, e.g. color, texture, interest points, motion (section 1.7).

6. Re-ranking. We post-process here the scores using the fact that videos statistically have an homogeneous content, at least locally (section 1.8).

This approach is quite similar to the one used by the IRIM group last year [START_REF] Delezoide | IRIM at TRECVID 2011: Semantic Indexing and Multimedia Instance Search[END_REF]. The main novelties are the inclusion of new descriptors, by some improvements in the pre-processins step and iprovements in the automatic fusion methods.

Descriptors

Thirteen IRIM participants (CEA-LIST, ETIS/LIP6, EURECOM, GIPSA, INRIA, LABRI, LIF, LIG, LSIS, LIRIS, Mines Télécom and Télécom ParisTech) provided a total of 128 descriptors, including variants of a same descriptors. These descriptors do not cover all types and variants but they include a significant number of different approaches including state of the art ones and more exploratory ones. The relative performance of these descriptors has been separately evaluated using a combination of LIG classifiers (see section 1.5). Here is a description of these descriptors (some of the descriptors used in 2011 are described in more details ih The TRECVid IRIM 2011 paper [START_REF] Delezoide | IRIM at TRECVID 2011: Semantic Indexing and Multimedia Instance Search[END_REF]):

CEALIST/tlep: texture local edge pattern [START_REF] Cheng | Image classification using color, texture and regions[END_REF] + color histogram 576 dimensions.

CEALIST/bov dsiftSC 8192: : bag of visterm [START_REF] Shabou | Locality-constrained and spatially regularized coding for scene categorization[END_REF]. Dense SIFT are extracted every 6 pixels. The codebook of size 1024 is built with K-means. The bag are generated with soft coding and max pooling. The final signature result from a three levels spatial pyramid 1024×(1+2×2+3×1) = 8192 dimensions: see [START_REF] Ballas | Aymen Shabou CEA LIST at TRECVID 2012: Semantic Indexing and Instance Search[END_REF] for details.

ETIS/global <feature>[<type>]x<size>:

(concatenated) histogram features [START_REF] Gosselin | Combining visual dictionary, kernel-based similarity and learning strategy for image category retrieval[END_REF], where: <type> can be:

<feature>
m1x1: histogram computed on the whole image m1x3: histogram for 3 vertical parts m2x2: histogram on 4 image parts <size> is the dictionary size, sometimes different from the final feature vector dimension.

For instance, with <type>=m1x3 and <size>=32, the final feature vector has 3 × 32 = 96 dimensions.

ETIS/vlat <desc type> dict<dict size> <size>: compact Vectors of Locally Aggregated Tensors (VLAT [START_REF] Negrel | Compact Tensor Based Image Representation for Similarity Search[END_REF]). <desc type> = low-level descriptors, for instance hog6s8 = dense histograms of gradient every 6 pixels, 88 pixels cells. <dict size> = size of the low-level descriptors dictionary. <size> = size of feature for one frame. Note: these features can be truncated. These features must be normalized to be efficient (e.g. L 2 unit length).

EUR/sm462: The Saliency Moments (SM) feature [START_REF] Redi | Saliency moments for image categorization[END_REF] is a holistic descriptor that embeds some locally-parsed information, namely the shape of the salient region, in a holistic representation of the scene, structurally similar to [START_REF] Oliva | Modeling the shape of the scene: A holistic representation of the spatial envelope[END_REF]. More details in [START_REF] Delezoide | IRIM at TRECVID 2011: Semantic Indexing and Multimedia Instance Search[END_REF]. [START_REF] Paris | Pyramidal Multi-level Features for the Robot Vision@ICPR 2010 Challenge[END_REF]. More details in [START_REF] Delezoide | IRIM at TRECVID 2011: Semantic Indexing and Multimedia Instance Search[END_REF].

INRIA
MTPT/superpixel color sift k1064: this visual feature is extracted based on superpixel segmentation [START_REF] Levinshtein | Turbopixels: Fast superpixels using geometric flows[END_REF]; it is an histogram combination of color and texture over superpixels of a given image; these histograms are computed basically using trained codebooks as bag-of-words.

Descriptor optimization

The descriptor optimization consists of two steps: power transformation and principal component analysis (PCA) reduction [START_REF] Safadi | Quaero at TRECVID 2011: Semantic Indexing and Collaborative Annotation[END_REF].

Power transformation: The goal of the power transformation is to normalize the distributions of the values, especially in the case of histogram components. It simply consists in applying an x ← x α (x ← -(-x) α if x < 0) tranformation on all components individually. The optimal value of α can be optimized by cross-validation and is often close to 0.5 for histogram-based descriptors.

Principal component analysis: The goal of PCA reduction is both to reduce the size (number of dimensions) of the descriptors and to improve performance by removing noisy components.

The optimization of the value of the α coefficient and of the number of components kept in the PCA reduction is optimized by two-fold cross-validation within the development set. In practice, it is done with the LIG KNNB classifier only (see section 1.4), since it is much faster when a large number of concepts (346 here) has to be considered and since it involves a large number of combinations to be evaluated. Trials with a restricted number of varied descriptors indicated that the optimal values for the kNN based classifier are close to the ones for the multi-SVM based one. Moreover, the overall performance is not very sensitive to the precise values for these hyper-parameters.

Classification

The LIG participant ran two types of classifiers on the contributed descriptors as well as their combination.

CEA LIST also run a classifier on its descriptors.

LIG KNNB:

The first classifier is kNN-based. It is directly designed for simultaneously classifying multiple concepts with a single nearest neighbor search. A score is computed for each concept and each test sample as a linear combinations of 1's for positive training samples and of 0's for negative training samples with weights chosen as a decreasing function of the distance between the test sample and the reference sample. As the nearest neighbor search is done only once for all concepts, this classifier is quite fast for the classification of a large number of concepts. It usually gives lower classification rates than the SVM-based one but is much faster.

LIG MSVM:

The second one is based on a multiple learner approach with SVMs. The multiple learner approach is well suited for the imbalanced data set problem [START_REF] Safadi | Evaluations of multi-learners approaches for concepts indexing in video documents[END_REF], which is the typical case in the TRECVID SIN task in which the ratio between the numbers of negative and positive training sample is generally higher than 100:1.

CEALIST LSVM:

The third one is a linear SVM classifier and was applied by CEALIST to its own high-dimensional descriptors.

LIG ALLC: Fusion between all available classifiers.

The fusion is simply done by averaging the classification scores produced by the two classifiers.

Their output is naturally (or by construction) normalized in the the [0:1] range. kNN computation is done using the KNNLSB package [START_REF] Quénot | KNNLSB: K Nearest Neighbors Linear Scan Baseline[END_REF]. Even though the LIG MSVM classifier is often significantly better than the LIG KNNB one, the fusion is most often even better, probably because they are very different in term of information type capture.

Evaluation of classifier-descriptors combinations

We evaluated a number of image descriptors for the indexing of the 346 TRECVID 2011 concepts. This has been done with two-fold cross-validation within the development set. We used the annotations provided by the TRECVID 2012 collaborative annotation organized by LIG and LIF [START_REF] Ayache | Video Corpus Annotation using Active Learning[END_REF]. The performance is measured by the inferred Mean Average Precision (MAP) computed on the 346 concepts. Results are presented for the two classifiers used, as well as for their fusion. Results are presented only for the best combinations of the descriptor optimization hyper-parameters.

Table 1 shows the two-fold cross-validation performance (trec eval MAP) within the development set and the performance (sample eval MAP) on the test set for most the descriptors (some variants are skipped) with the LIG ALLC classifier combination; dim is the original number of dimensions of the descriptor vector, exp is the optimal value of the α coefficient, Pdim is the number of dimensions of the descriptor vector kept after PCA reduction.

At the time of the writing, not all descriptor × classifier combinations were computed. For some descriptors, only the scores from the KNN classifier was available and are displayed in the table (KNN only) and/or som are not available on the test set. Some more were not available at the time of run submissions and not all of them were included in the fusion or some less performing versions were used.

Performance improvement by fusion of descriptor variants and classifier variants

As in previous years, we started by fusing classification scores from different variants of a same descriptor and from different classifiers of a same variant of a same descriptor. This is done as first levels of hierarchical late fusion, the last ones being done using dedicated methods as described in section 1.7. Three levels are considered when applicable: fusions of different classifiers of a same variant of a same descriptor, fusion of different variants of a same descriptor according to a dictionary size, and fusion of different variants of a same descriptor according to a pyramidal decomposition. While the last levels of fusion attempt to improve the overall performance by fusing information of different types (e.g;. color, texture, percepts or SIFT), the first fusion levels attempt to improve the robustness of the classification from a given type. More details on this approach can be found in the previous TRECVid IRIM papers [START_REF] Gorisse | IRIM at TRECVID 2010: High Level Feature Extraction and Instance Search[END_REF][START_REF] Delezoide | IRIM at TRECVID 2011: Semantic Indexing and Multimedia Instance Search[END_REF].

Final fusion

Two IRIM participants (LISTIC and LIMSI) worked on the automatic fusion of the classification results.

The fusion started with the original classification scores and/or with the results of previous fusions of descriptor variants and/or classifier variants as described in the previous section. A comparison of the LISTIC and LIMSI automatic fusion methods, along with another fusion method tried in the context of the Quaero group using some of the same classification results, and an arithmetic mean and the best attribute per concept, is given in [START_REF] Strat | Hierarchical Late Fusion for Concept Detection in Videos[END_REF].

LISTIC fusion

We will call an "attribute" the set of scores obtained by applying any of the KNNB, MSVM or ALLC classifiers on a descriptor. The inputs for this fusion approach are the KNNB, MSVM and ALLC scores obtained from various descriptors. We normally use the ALLC scores, but if they are not available for a particular descriptor, the KNNB and/or MSVM scores for that descriptor are used instead.

As explained in [START_REF] Strat | Hierarchical Late Fusion for Concept Detection in Videos[END_REF], this fusion approach treats each semantic concept independently. It automatically filters out irrelevant attributes for that concept, then it automatically groups highly-correlated attributes in an iterative manner. The method does not take into account the type of descriptor or the type of classifier for generating groupings, instead it groups based on the correlation of output scores. Being an automatic fusion method, it can easily be extended to large sets of available attributes, without having to manually specify groupings. The method consists in the following steps:

1. determine the individual relevance of each attribute for the target concept. The relevance is taken as the average precision α for that concept on the training dataset, normalized by the proportion of true positives in the training dataset.

2. retain only attributes with a relevance higher than 1 (better than random classification). Additionally, the attributes must have at least 1/8th of the relevance of the best one, so as not to "pollute" the good attributes with bad ones.

3. Some of the retained attributes are highly correlated, so we look for the pair with the maximum correlation and fuse it into a single attribute through arithmetic mean. We keep track of how many original attributes have already been fused into each of the members of the pair, in order to balance the members of the pair correctly. After replacing the pair with the new attribute, we update the correlation between the resulting attribute and the remaining ones.

4. The previous step is repeated many times, until a sufficiently correlated pair can no longer be found. This has a dimensionality-reduction effect and also helps to reduce the classification "noise".

The correlation measure used is the correlation coefficient of the raw (unnormalized) classification scores. We consider a pair of attributes as correlated if (a) the correlation coefficient for all video shots is at least 0.75, to ensure that the two attributes give similar information on a global scale, and (b), the correlation between the scores for just the positive shots must be at least 0.65, to ensure that the positives tend to be classified in the same way. We add the second constraint because in TRECVID, most of the target concepts have very few positives, and otherwise, the classification scores on the negatives would dominate the correlation. Now, the resulting attributes are again filtered based on their average precision on the training set, using the same criteria as before. Afterwards, the individual relevances of each remaining attribute are used as weights for a weighted arithmetic mean, thus obtaining the final classification score. In the end, because all the previous steps consisted of selections and averaging, without any normalization operations, this approach is in fact a weighted arithmetic mean, with a more elaborate way of choosing the weights.

This method has given the TRECVid Semantic Indexing Full Task runs F A IRIM1 1 (which also uses a temporal re-ranking of video shots) and F A IRIM3 3 (without temporal re-ranking).

LIMSI community-driven hierarchical fusion

Let K be the number of available classifiers and N the number of video shots. Each classifier k ∈ {1 . . . K} provides scores x k = [x k1 , . . . , x kN ] indicating the likelihood for each shot n ∈ {1 . . . N} to contain the requested concept. The objective is to find a combination function f so that the resulting classifier x = f (x 1 , . . . , x K ) is better than any of its components, and as good as possible.

Graph of classifiers

Let us denote ρ ij the Spearman rank correlation coefficient of two classifiers i and j. We then define the agreement A ij between two classifiers i and j as A ij = max (0, ρ ij ).

A complete undirected graph G is constructed with one node per classifier. Each pair of classifiers (i, j) is connected by an undirected edge, whose weight is directly proportional to A ij . Based on this graph G, classifiers can be automatically grouped into communities using the so-called Louvain approach proposed by Blondel et al. [START_REF] Blondel | Fast Unfolding of Community Hierarchies in Large Networks[END_REF]. Step 1: community detection. Classifiers are automatically grouped into C communities using the Louvain method described above;

Step 2: intra-community fusion. Classifiers from each community are combined by simple sum of normalized scores, in order to obtain one new classifier per community (classifiers A to E in Figure 1):

x c = k=K k=1 δ c (k) x k with δ c (k) = 1 if classifier k is part
of community c (and 0 otherwise);

Step 3: inter-community fusion. Those new classifiers are then combined using weighted sum fusion of normalized scores: x = c=C c=1 α c x c . To this end, the performance α c (average precision) of each of these new community classifiers needs to be estimated using a development set.

F A IRIM2 2: re-ranked version of F A IRIM4 4 using the method described in section 1.8.

Re-ranking

Video retrieval can be done by ranking the samples according to their probability scores that were predicted by classifiers. It is often possible to improve the retrieval performance by re-ranking the samples. Safadi and Quénot in [START_REF] Safadi | Re-ranking by Local Re-scoring for Video Indexing and Retrieval[END_REF] propose a re-ranking method that improves the performance of semantic video indexing and retrieval, by re-evaluating the scores of the shots by the homogeneity and the nature of the video they belong to. Compared to previous works, the proposed method provides a framework for the re-ranking via the homogeneous distribution of video shots content in a temporal sequence. The experimental results showed that the proposed re-ranking method was able to improve the system performance by about 18% in average on the TRECVID 2010 semantic indexing task, videos collection with homogeneous contents. For TRECVID 2008, in the case of collections of videos with non-homogeneous contents, the system performance was improved by about 11-13%.

Evaluation of the submitted runs

IRIM officially submitted the four F A IRIM1 1 to F A IRIM4 4 runs that are described in section 1.7. Table 2 presents the result obtained by the four runs submitted as well as the best and media runs for comparison. The best IRIM run corresponds to a rank of 4 within the 16 participants to the TRECVID 2012 full SIN task. [START_REF] Over | TRECVID 2012 -An Overview of the Goals, Tasks, Data, Eval[END_REF] consists in finding segments of videos in the data set which contain instances of these entities. Each instance is represented by a few example images. Hence if we can consider the set of video clips as a visual database, the problem consists in retrieval of each instance in this database.

Global approach

To represent the clips we extract several keyframes of each individual video clip. For a given instance, we use each example image, from the available set, as a query image. We compute a similarity between this query example image and the keyframes of all video clips. We then produce a intermediary result where we have the similarity S e,i,k,c between each example image (e) of each instance (i) and each keyframe (k) of each video clip (c). We then have to fuse these intermediary results to obtain a final result that is similarity S i,c between each instance (i) and each clip (c) Within the IRIM consortium, several methods of four members (CEA, CNAM, LaBRI, LISTIC) were tested and their results fused.

Members methods

CEA Markrs

The Markrs methods follows the wellknown framework of keypoint matching described in [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF][START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. Here, the SURF scheme [START_REF] Bay | SURF: Speeded Up Robust Features[END_REF] was used for keypoints detection and descriptors computations. Then SURF description is quantized into integer values in [0, 255], leading to a compact description for each keypoint in less than 80 bytes.

Two filtering steps are added to reduce bad keypoints matches. First, matches must pass the test of relative nearest-neighbors proposed by D. Lowe in [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF]. The second step selects within the previous results those that provides a similar geometrical configuration of keypoints in the query-candidate couple of images. We consider only simple similarities, and not complete homographies, that are faster to compute. The final result list is composed of images having more than p keypoints fitting the geometrical model (p ≥ 5).

CEA Snow

The Snow method is based on a wellknown color histogram for global image description. The color histogram counts the occurences of 162 shades in the HSV color space. The similarity between two images is measured as the inverse of a dLog distance between two histograms, defined as: [START_REF] Stehling | A compact and efficient image retrieval approach based on border/interior pixel classification In 11[END_REF]:

dLog(q, d) = i<M i=0 |f (q[i]) -f (d[i])| (1)
f (x) =    0, if x = 0 < α 1, if 0 < x ≤ 1 ⌈log 2 x⌉ + 1 otherwise (2)
Where q and d are two histograms with M bins and ⌈.⌉ is the ceilling function (). These methods were also used in individual CEA LIST submission [START_REF] Ballas | Aymen Shabou CEA LIST at TRECVID 2012: Semantic Indexing and Instance Search[END_REF].

CNAM This method uses inverted lists for pairs of visual words. The frames of the videos in the database are described with local features, then visual words are obtained using a vocabulary tree and an inverted list is created for every pair of visual words co-occurring in at least one frame. The inverted lists include geometric information relating the scales of the two features. Every query is described in the same way and retrieval is performed by accessing the selected inverted lists. This solution was inspired by the method put forward in [START_REF] Poullot | Indexing local configurations of features for scalable content-based video copy detection[END_REF] that used triplets of Harris features (and associated geometric information) for fast content-based video copy detection. Since the TRECVID INS data shows scale and viewpoint variations, we employed instead SURF features. Individual SURF features include data regarding scale and orientation, which allows to use pairs of features rather than triplets, with an expected positive impact on recall, while keeping geometric information to improve prcision. Retrieved frames are ranked according to a similarity measure that takes into account the number of pairs from the query that are present in the frame, the similarity of corresponding individual SURF features (based on the path in the vocabulary tree), the ratios keypoints scale to distance between keypoints, and the relative orientations of the keypoints. The use of pairs of keypoints together with simple geometric information supports other possibilities that were not fully explored in the current implementation. Retrieval is indeed very fast: even when the database is composed of all the frames (and not only of the keyframes), less than 15 minutes are required for obtaining all the results to a query.

LaBRI We used the baseline Bag-Of-Visual-Words (BOVW) model based on interest point descriptor [START_REF] Sivic | Video google: a text retrieval approach to object matching in videos[END_REF]. It is used both for object and for frame signature construction. The descriptor used is SURF [START_REF] Bay | SURF: Speeded Up Robust Features[END_REF]. The unsupervised clustering K-means++ with a large number of clusters (16K), with the L2 distance, was used for dictionary computation. The complement of histogram intersection was used to compare signatures. These methods will be referenced as LaBRI Bow Obj and LaBRI Bow Frm.

For frame-based queries, we also used BOVW with affine deformation of object mask, described in [START_REF] Mansencal | Domenger Search of objects of interest in videos[END_REF]. This method uses a correlation kernel, deforming object mask according to Pan/Tilt/Zoom affine model. The correlation was done by full search in the affine parameter space. Pan and Tilt parameters were chosen in such a way that query instance mask overlaps the DB frame at least two third of its area. The Zoom factor were chosen from the set 0.25, 0.5, 1, 2. This method is obviously more computationally demanding than the traditional BOVW. Indeed, signatures can not be computed in a processing step for all the images of the DB, but have to be computed in image area overlapped by image mask. This method will be referenced as LaBRI Bow Obj Aff.

A second approach is a Bag-Of-Regions (BOR) model, as proposed by Vieux et al. in [START_REF] Vieux | Content based image retrieval using bag of region[END_REF], that extends the traditional notion of BOVW vocabulary to region based [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF]. The HSV histogram was computed as region feature. We set a uniform quantizing parameters in order to limit the feature size to approximately 100 bins (45+32+32) and to privilege the finest encoding of Hue component. We used the incremental clustering algorithm described in [START_REF] Lughofer | Extensions of vector quantization for incremental clustering[END_REF] and modified in [START_REF] Vieux | Content based image retrieval using bag of region[END_REF], with 2K clusters and L2 distance. Region-based approach was deployed only for frame signatures construction. To compare signatures, we used the L1 distance. This method will be referenced as LaBRI Bor Frm.

For development and test sets, we computed their proper codebooks as we are not granted that the two sets have the same distribution in proposed description spaces.

LISTIC The two methods LISTIC SIFT L2 k and LISTIC SIFT L2 retina k used in the SIN task were also used in the INS task, with k=1024 and complement of histogram intersection for signatures comparison. Codebook computed on dataset of SIN task was used.

Fusion

Each described members method was used to produce intermediary results. Thus for each method (m) , we have a similarity S m,e,i,k,c between each example image (e) of each instance (i) and each keyframe (k) of each video clip (c). We have to fuse these similarities to obtain a similarity for an instance (i) and a clip (c).

We used a limited number combination operators: CombMAX, CombSUM, CombANZ, CombMNZ [START_REF] Fox | Combination of Multiple searches[END_REF], defined in table 3.

We have tested two late fusion schemes. A truly late fusion scheme considers all the similarities S m,e,i,k,c at once. In a two-step late fusion scheme, we first merge the results for a given method (m), and then globally. Besides, weights can be used to give an asymmetric importance to the various intermediary results. Here, all intermediary results have been previously normalized. These two fusion schemes are described by the equations 3 and 4, where α m and β m are weights that sum to 1.

S i,c = Comb 1 (α m * S m,e,i,k,c ) (3)
S i,c = Comb 2 (β m * S e,i,k,c ) with S e,i,k,c = Comb m (α m * S m,e,i,k,c ) (4)
A Combination operator will be noted Comb[S] if applied to score, and Comb[R] if applied to rank. We have tested several combination operators with these two fusion schemes, applied both to score and to rank. We have also tested with a limited combination of weights. The best results were obtained with the two-step fusion scheme. Moreover, as performance of various methods is not homogeneous, we tried to find the best combination operator and the similarity to use for each individual method, both for 2010 and 2011 queries and datasets. Theses choices are presented in table 4.

We tested a limited number of combinations of weights. Finally, we used three decreasing functions. For N combined intermediary results, we used the formulae in ( 5), ( 6), [START_REF] Oliva | Modeling the shape of the scene: A holistic representation of the spatial envelope[END_REF] for the i-th intermediary result, ranked from left to right, within the fusion step.

β1(i) = 1 -i/N (5) β2(i) = 1/(1.4 i ) (6) β3(i) = 1/(3 i ) (7)
Various experiments on 2010 and 2011 datasets showed that CEA marks indivudually was giving the best results compared to all other methods. BOVW and BOR methods gave inferior, but promising results on 2011 dataset. So these results were given a bigger weight in the second fusion step, in eq. ( 4).

If we consider all possibilities of fusion: two-step or one-step ((3), ( 4)), combination operators (3), weight functions (eq. ( 5), ( 6), ( 7)), choice of rank (R) or similarity (S), and order the results of individual methods according to decrease of performances, we get 1536 possibilities for one run. Obviously, it was impossible to test all these combinations. Hence, we proposed 4 runs on the basis of a two-step fusion scheme, all three weight functions and a best fusion operator for individual methods. The choice of rank (R) or similarity (S) for the final fusion step was done on the basis of tests on 2010 and 2011 sets.

The submitted runs were:

Run1 = CombM AX[R]( β1(0) * CombSU M [S](CEA markrs), β1 (

Discussion

MAP and rank of our runs are presented in The results are very much close. Nevertheless, we can see that an appropriate fusion odes bring an improvement (run IRIM 2 vs CeaList 2). In the future while improving individual methods, we plan also to explore different fusion operators , such as described in [START_REF] Csurka | An empirical study of fusion operators for multimodal image retrieval[END_REF].

Data sharing

We propose to reuse and extend the organization that has been developed over four years within the members of the IRIM project of the French ISIS national Research Group (see [START_REF] Delezoide | IRIM at TRECVID 2011: Semantic Indexing and Multimedia Instance Search[END_REF] and section 1 of this paper).

It is based on a limited number of simple data formats and on a (quite) simple directory organization. It also comes with a few scripts and procedures as well as with some sections for reporting intermediate results. The supporting structure is composed of a wiki (http://mrim.imag.fr/trecvid/wiki) and a data repository (http://mrim.imag.fr/trecvid/sin12). The wiki can be accessed using the TRECVid 2012 active participant username and password and the data repository can be accessed using the TRECVid 2012 IACC collection username and password.

A general rule about the sharing of elements is that:

• any group can share any element he think could be useful to others with possibly an associated citation of a paper describing how it was produced;

• any group can use any element shared by any other group provided that this other group is properly cited in any paper presenting results obtained using the considered element, exactly as this was the case in the previous years for the shared elements like shot segmentation, ASR transcript or collaborative annotation. Groups sharing elements get "rewarded" via citations when their elements are used.

Shared elements can be for instance: shot or key frame descriptors, classification results, fusion results. For initiating the process, most IRIM participants agreed to share their descriptors. Most classification and fusion results obtained are also shared. These should also be available on the 2013 TRECVID SIN collection. Descriptors, classification scores or fusion results from other TRECVid particpants are most welcome. See the wiki for how to proceed.

  is chosen among lab and qw: lab: CIE L*a*b* colors qw: quaternionic wavelets (3 scales, 3 orientations)

Figure 1 :

 1 Figure 1: Community-driven hierarchical fusion F A IRIM4 4: Hierarchical fusion It can be divided into three consecutive steps:

1 ) 8 ) 9 )

 189 * CombAN Z[S](LaBRI Bow Obj), β1(2) * CombM AX[S](LaBRI Bow Obj Af f ), β1(3) * CombM AX[S](LIST IC SIF T L2 1024), β1(4) * CombAN Z[S](LaBRI Bow F rm), β1(5) * CombM AX[R](CEA snow), β1(6) * CombAN Z[S](LaBRI Bor F rm), β1(7) * CombAN Z[S](LIST IC SIF T L2 retina 1024), β1(8) * CombM AX[S](CN AM )) (Run2 = CombSU M [R]( β3(0) * CombSU M [S](CEA markrs), β3(1) * CombAN Z[S](LaBRI Bow Obj), β3(2) * CombM AX[S](LaBRI Bow Obj Af f ), β3(3) * CombM AX[S](LIST IC SIF T L2 1024), β3(4) * CombAN Z[S](LaBRI Bow F rm), β3(5) * CombM AX[R](CEA snow), β3(6) * CombAN Z[S](LaBRI Bor F rm), β3(7) * CombAN Z[S](LIST IC SIF T L2 retina 1024), β3(8) * CombM AX[S](CN AM )) (Run3 = CombM AX[R]( β2(0) * CombSU M [S](CEA markrs), β2(1)* CombM AX[S](LaBRI Bow Obj Af f ), β2(2) * CombAN Z[S](LaBRI Bor F rm), β2(3)* CombM AX[R](CEA snow))(10)Run4 = CombM AX[R]( β2(0) * CombSU M [S](CEA markrs), β2(1) * CombAN Z[S](LaBRI Bow Obj), β2(2) * CombM AX[S](LaBRI Bow Obj Af f ), β2(3) * CombAN Z[S](LaBRI Bow F rm), β2(4) * CombM AX[R](CEA snow), β2(5) * CombAN Z[S](LaBRI Bor F rm))

Table 1 :

 1 Performance of the classifier and descriptor combinations

	Descriptor

Table 2 :

 2 InfMAP result and rank on the test set for all the 46 TRECVID 2012 evaluated concepts (full task).

	System/run	MAP rank
	Best run	0.3210	1
	F A IRIM1 1 0.2379	12
	F A IRIM3 3 0.2278	13
	F A IRIM2 2 0.2248	16
	F A IRIM4 4 0.2153	18
	Median run	0.1944	26
	2 Instance Search
	Given visual examples of entities of limited number of
	types: person, character, object or location, Instance
	Search (INS) task

Table 3 :

 3 Definitions of different combination operators descriptors. Regions in image plane are obtained by segmenting images by Felzenszwalb and Huttenlocher method

Table 4 :

 4 Best combination operator with similarity type used for each individual methods

table 5 .

 5 The two-step combination method applied to Markrs

	Run	MAP	Rank/79
	IRIM 2	0.1192	29
	IRIM 4	0.1173	31
	IRIM 1	0.1171	32
	IRIM 3	0.1162	33
	CeaList 2 0.1135	34

Table 5 :

 5 MAP and rank of 4 IRIM runs and CeaList 2 runs on test set 2012, among 79 fully automatic submitted runs method was also submitted by CEA-LIST as an individual run, as CeaList 2. It is added for reference.

Acknowledgments

This work has been carried out in the context of the IRIM (Indexation et Recherche d'Information Multimédia) of the GDR-ISIS research network from CNRS. Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see https://www.grid5000.fr).