Denys Dutykh
email: denys.dutykh@univ-savoie.fr

Efim Pelinovsky
email: pelinovsky@hydro.appl.sci-nnov.ru

NUMERICAL SIMULATION OF A SOLITONIC GAS IN KDV AND KDV-BBM EQUATIONS

Keywords: solitonic gas, KdV equation, BBM equation, statistical description, statistical moments, Stokes-Ursell number 2010 Mathematics Subject Classification. 76B25, 76B15, 35Q53

published or not. The documents may come

Numerical simulation of a solitonic gas in KdV and

KdV-BBM equations Denys Dutykh, Efim

Introduction

Solitary wave solutions play the central rôle in various nonlinear sciences ranging from hydrodynamics to solid and plasma physics [44,[START_REF] Salupere | On the long-time behaviour of soliton ensembles[END_REF]30]. These solutions can propagate without changing its shape. However, the most intriguing part consists in how these solutions interact with each other. The binary interactions of solitary waves have been studied in the context of various nonlinear wave equations [44,23,[START_REF] Craig | Solitary water wave interactions[END_REF]25,9,6]. It is well known that in integrable models the collision of two solitons is elastic, i.e. they interact without emitting any radiation. In non-integrable models usually the interactions are nearly elastic [4].

The collective behaviour of soliton ensembles is much less understood nowadays. When a large number of solitary waves are considered simultaneously the researchers usually speak about the so-called solitonic turbulence or a solitonic gas. The literature on this topic is abundant. Some recent studies on solitonic gas turbulence in the KdV framework include [START_REF] Pelinovsky | Numerical modeling of the KdV random wave field[END_REF]35,[START_REF] Pelinovsky | Two-soliton interaction as an elementary act of soliton turbulence in integrable systems[END_REF]. The solitonic turbulence in nonintegrable NLS-type equations was studied in [START_REF] Zakharov | Soliton turbulence[END_REF][START_REF] Dyachenko | Soliton turbulence in nonintegrable wave systems[END_REF] and the authors showed that in conservative nonintegrable systems the solitonic gas is a statistical attractor whose dimension decreases with time. Recently it was shown both numerically and experimentally that solitonic ensembles appear in the laminar-turbulent transition in a fibre laser [START_REF] Turitsyna | The laminar-turbulent transition in a fibre laser[END_REF], modeled by a non-integrable nonlinear Schrödinger-type equation. However, the dominant number of studies is based on integrable models. This apparent contradiction motivated mainly our investigation to quantify the non-integrable effects onto the collective behaviour of solitons.

An approximate theoretical description of solitonic gases was proposed by V. [START_REF] Zakharov | Kinetic Equation for Solitons[END_REF] [45] using the kinetic theory. Later this research direction has been successfully pursued by G. El & A. Kamchatnov [START_REF] El | Kinetic Equation for a Dense Soliton Gas[END_REF][START_REF] El | Kinetic Equation for a Soliton Gas and Its Hydrodynamic Reductions[END_REF] who used the Inverse Scattering Technique (IST) [START_REF] Ablowitz | Solitons and the Inverse Scattering Transform[END_REF] limited only to the integrable models. In this study the problem of solitonic gases will be investigated using the methods of direct numerical simulation. The evolution of random wave fields including solitonic gases was simulated numerically in [START_REF] Pelinovsky | Numerical modeling of the KdV random wave field[END_REF]35,[START_REF] Dutykh | Geometric numerical schemes for the KdV equation[END_REF] using symplectic, multi-symplectic and pseudo-spectral methods. However, previous investigators considered only a limited number of solitons (a few dozens) to simulate a solitonic gas. In this study we will adopt the pseudo-spectral method since it provides the high accuracy and computational efficiency necessary to handle large computational domains. Our goal will consist in:

• investigate the influence of soliton interactions on statistical characteristics of the wave field, • construct the Probability Density Function (PDF) and compute the first four statistical moments of the solitonic turbulence, • study the rôle of non-integrable terms on the characteristics of soliton ensembles.

The present manuscript is organized as follows. In Section 2 we derive the governing equation used in this study and in Section 3 numerical results on a solitonic gas dynamics are presented. Finally, the main conclusions of this study are outlined in Section 4.

Mathematical model

As the starting point we choose the celebrated Korteweg-de Vries equation [22,24,[START_REF] Johnson | Camassa-Holm, Korteweg-de Vries and related models for water waves[END_REF][START_REF] Pelinovsky | Numerical modeling of the KdV random wave field[END_REF] (in dimensional variables) which models the undirectional propagation (here in the rightwards direction) of weakly nonlinear and weakly dispersive waves:

η t + c 1 + 3h 2 η η x + ch 2 6 η xxx = 0, (2.1)
where η(x, t) is the vertical excursion of the free surface above the still water level, h is the uniform undisturbed water depth and c = √ gh is the speed of linear gravity waves (g being the gravity acceleration).

The KdV equation (2.2) is known to be integrable [18,[START_REF] Chan | A computer study of finite-amplitude water waves[END_REF]. However, the full water wave problem is known to be a non-integrable system, since the interaction of solitary waves is inelastic [7, [START_REF] Fenton | A Fourier method for solving nonlinear water-wave problems: application to solitary-wave interactions[END_REF]8,9]. Moreover, we will modify the original equation (2.2) in order to include an additional dispersive term of the BBM-type [START_REF] Benjamin | Model equations for long waves in nonlinear dispersive systems[END_REF]. Thus, the resulting KdV-BBM equation will not be integrable [START_REF] Francius | Wave dynamics in nonlinear media with two dispersionless limits for long and short waves[END_REF][START_REF] Dutykh | Finite volume methods for unidirectional dispersive wave models[END_REF]21].

Consider the following scaled dependent and independent variables:

η ← η a 0 , x ← x l , t ← ct l ,
where a 0 is the characteristic wave amplitude and l is the characteristic wavelength. In dimensionless variables KdV equation (2.2) reads:

η t + 1 + 3ε 2 η η x + µ 2 6 η xxx = 0,
where parameter ε := a 0 /h measures the nonlinearity and µ 2 := h /l 2 is the dispersion parameter. The relative importance of these two effects is measured by the so-called Stokes-Ursell number [42] (sometime denoted as Ur):

S := ε µ 2 ≡ a 0 l 2 h 3 .
The last equation can be further simplified if we perform an additional change of variables:

η ← 3µ 2 2 η, x ← √ 6 µ (x -t), t ← √ 6 µ t,
which yields the following simple equation including explicitly the Stokes-Ursell number S:

η t + Sηη x + η xxx = 0. (2.
2) The last scaled KdV equation can be further generalized by using the low-order asymptotic relations in order to alternate higher-order terms as it was proposed by [START_REF] Bona | A model for the two-way propagation of water waves in a channel[END_REF] [5] and Nwogu (1993) [28]. This step is rather standard and we do not provide here the details of the derivation:

η t + Sηη x + η xxx -δη xxt = 0, (2.3)
where δ ∈ R is a free parameter. . Below we will study the solitonic gas behaviour1 under the dynamics of the KdV (2.2) (δ = 0) and KdV-BBM (2.3) (δ = 0) equations. In the absence of the KdV term, we recover the celebrated BBM equation derived in [START_REF] Peregrine | Calculations of the development of an undular bore[END_REF][START_REF] Benjamin | Model equations for long waves in nonlinear dispersive systems[END_REF].

Remark 1. We note that for a particular value of the Stokes-Ursell number S ≡ 1 another simpler scaling is possible when all the lengthes (x and η) are scaled by the mean water depth h. However, we do not adopt it in this study since below in Section 3.1 the dependence of some statistical characteristics on the Stokes-Ursell number S is investigated.

Properties 2.1.1. Linear well-posedness

In order to ensure the linear well-posedness of equation ( 2.3), the free parameter δ has to satisfy the following constraint δ ≥ 0. In the following we will consider only nonnegative values of this parameter. Recall that for δ = 0 we recover an scaled version of the classical KdV equation (2.2).

Ivariants

Assuming that the solution η(x, t) has either the compact support or decays sufficiently fast at the infinity along with its first derivative (η → 0, η x → 0 as x → ±∞), one can easily show that the following quantities are conserved [START_REF] Dutykh | Finite volume methods for unidirectional dispersive wave models[END_REF]:

I 1 (t) = ˆR η(x, t) dx, I 2 (t) = ˆR η 2 (x, t) + δη 2 x (x, t) dx.
In other words, I 1 (t) ≡ I 1 (0) and I 2 (t) ≡ I 2 (0), ∀t > 0. The invariant I 1 (t) is related to the mass conservation property, while the integral I 2 (t) can be assimilated to the generalized kinetic energy. The conservation of these quantities has not only the theoretical importance, but also the practical one. For example, it will allow us to check the global accuracy of the employed numerical scheme. We note also that the same invariants hold also in finite, but periodic domains (below we use periodic boundary conditions). The conservation laws to the BBM equation can be found in [29].

Solitary wave solutions

Equation (2.3) admits an exact localized (solitary) travelling wave solution which can be found analytically:

η(x, t) = a sech 2 1 2 κ(x -c s t) , κ := aS 3 + aSδ , c s := 1 3 aS. (2.4)
The dependence of the solitary wave shape on parameters a, S and δ is shown on Figure 1.

For instance, one can see that solitary waves become thiner when the amplitude a (and hence the speed) and/or the Stokes-Ursell number S are increased (see Figure 1(a,c)). On the other hand, the increase of the BBM coefficient δ leads to the growth of the tail (see Figure 1(b)). If both parameters S and δ are increased simultaneously, the 'thinning' effect of the Stokes-Ursell number dominates (see Figure 1(d)). The solitary waves interact elastically only in the integrable KdV case (δ ≡ 0) [26], while some dispersive tails are generated after the interaction in the general KdV-BBM model. We note that some authors came to the wrong conclusion about the elasticity of solitary wave interactions in the BBM equation [START_REF] Eilbeck | Numerical study of the regularized long-wave equation. II: Interaction of solitary waves[END_REF] on the basis of low accuracy numerical simulations.

Numerical results

In order to solve numerically equation (2.3) we use a Fourier-type pseudo-spectral method with 3/2-antialiasing rule [START_REF] Trefethen | Spectral methods in MatLab[END_REF]. embedded adaptive 9(8) Runge-Kutta scheme [43]. The time step is chosen adaptively using the so-called H211b digital filter [START_REF] Söderlind | Digital filters in adaptive time-stepping[END_REF][START_REF] Cooker | Reflection of a high-amplitude solitary wave at a vertical wall[END_REF] to meet some prescribed error tolerance (generally of the order of machine precision ∼ 10 -15 ). The number of Fourier modes, the length of the computational domain and other numerical parameters are specified in Table 1.

In long time simulations presented below the invariants I 1,2 were conserved within the numerical accuracy 10 -11 and 10 -9 correspondingly. For the sake of illustration the error of the invariant I 1 (t) conservation in KdV and KdV-BBM simulations is shown on Figure 2. This accuracy is satisfactory to draw some robust conclusions on a solitonic gas statistical behaviour.

The initial condition for the KdV equation is composed of a finite number N s of solitons with random ampmlitudes a i ∼ N(a, σ), i = 1, . . . , N s separated by quasi-uniform distance ∆x s which is randomized to improve the ergodocity of the initial state:

∆x s := ∆x s + N(0, σ 2 ),
where ∆x s denotes the mean value reported in Table 1 and N(µ, σ) is the normal distribution with the mean µ and the variance σ. The solitonic gas initial state generated in this way is depicted on Figure 3(a). We use the same parameters for the initial solitonic gas state for the simulation with the KdV-BBM equation except for the domain size and the spacing between two solitons (see Table 1). They are larger in the KdV-BBM case since the soliton width increases with the parameter δ ≥ 0 (see Figure 1). Consequently, we can say that two initial conditions are approximatively isomorphic up to the horizontal coordinate stretching transformation. The simulation times T 1 (KdV) and T 2 (Kdv-BBM) are chosen so that an average soliton has enough time to go around the whole computational domain. The final states of the simulations are shown on Figures 3(b,c). One can notice how the initially quasi-uniform distribution of solitons is mixed by forming instantaneous soliton clusters as long as some void spaces due to the mass conservation property. The complete space-time dynamics simulated using the KdV and KdV-BBM equations is depicted on Figures 4(a,b). Individual lines correspond to solitons trajectories. The convergence of these lines corresponds to solitons collisions. It might appear on Figure 4 that collisions involve multiple solitons, however it is not the case. A zoom on a portion of the space-time domain is shown on Figure 5. One can see that interactions are only binary in agreement with [45]. It is interesting to estimate the total number of collisions in our simulation. The models under consideration are unidirectional. Thus, we can experience only overtaking collisions. Since the simulation time is chosen so that almost every soliton has enough time to travel across the whole computational domain, the number of collisions scales with O(m + m -), where m ± is the number of solitons which travel with the speed above (below) the average. By construction of the initial condition we have m ± = O( 1 2 N s ). Consequently, the total number of collisions scales with O( 1 4 N 2 s ). Since the initial conditions are approximatively self-similar, after an appropriate rescaling of the spatial and time variables, the space-time dynamics is similar in both simulations (see Figures 4(a,b)). The difference between two simulations can be noticed if one makes a zoom on solitons background in order to see small radiating oscillations due to the inelasticity of collisions in the KdV-BBM case. This zoom on a portion of the computational domain is shown on Figure 6. In contrast, Figure 6(a) shows the absence of phonon modes in the KdV simulation.

It is custom to use the statistical methods to describe random wave fields [START_REF] Boccotti | Wave Mechanics for Ocean Engineering[END_REF]19]. The probability distribution of the normalized free surface elevation η 0 (x, t) := η(x, t)η / η 2 1/2 at times t = 0 and t = T 1 is shown on Figure 7 (under the KdV dynamics). One can see that this distribution is quasi-stationary which is a direct consequence of the fact that solitons preserve perfectly their shape during the interactions. We note that the KdV-BBM numerical result shows the same invariance property since the inelasticity Please, note that the final simulation times T 1 = T 2 (see Table 1 for the values of T 1,2 ). Parameter τ 0 denotes the time needed for an average soltion to go over the whole computational domain.

Estimation of statistical moments

In order to estimate efficiently kurtosis and skewness for various values of parameters S and δ without performing a series of direct numerical simulations we will adopt an approximate analytical method of statistical moments estimation employed also in [36].

Let us introduce the average density of a solitonic gas:

ρ := N s 2ℓ ,
where N S is the number of solitons and ℓ is the half-length of the computational domain (see also Table 1). We will assume that the solitonic gas is rarefied, i.e. ρ ≪ 1. Under this assumption we can represent approximatively the instantaneous free surface elevation η(x, t) as a linear superposition of distinct solitary waves (the interacting part is neglected):

η(x, t) ≈ Ns i=1 η i (x, t) = Ns i=1 a i sech 2 1 2 κ i ξ , ξ := x -c si t -x i ,
where {a i }, {c si } and {x i } are respectively the amplitudes, speeds and phase shifts of individual solitary waves. By assuming that the supports of solitons do not overlap, we can estimate the statistical moments of any order. For the sake of simplicity let us consider the first order µ 1 , i.e. the mean:

µ 1 = η(x, t) = 1 2ℓ Ns i=1 ˆℓ/2 -ℓ/2 η i (ξ) dξ ≈ 1 2ℓ Ns i=1 a i ˆ+∞ -∞ sech 2 1 2 κ i ξ dξ = 4ρ a i κ i ,
where we took the limit ℓ → +∞ in order to compute analytically the integrals. Note also that the solitary wave parameter κ i is a function of the amplitude a i according to the relations given in (2.4). Higher order moments µ n = η n , n > 1 can be computed in a similar way. In this study we will need the moments up to the fourth order: Using these moments we can estimate the skewness ς and kurtosis κ in the rarefied gas limit ρ → 0:

µ 2 = 8 3 ρ a 2 i κ i , µ 3 = 32 15 ρ a 3 i κ i , µ 4 = 64 35 ρ a 4 i κ i . 0 0.
ς = (η -µ 1 ) 3 (η -µ 1 ) 2 3/2 = µ 3 µ 3/2 2 + O(ρ 1/2 ) ≈ 2 √ 3 5 √ 2 ρ -1/2 a 3 i /κ i a 2 i /κ i 3/2 , (3.1)
κ = (η -µ 1 ) 4 (η -µ 1 ) 2 2 = µ 4 µ 2 2 + O(1) ≈ 9 35 ρ -1 a 4 i /κ i a 2 i /κ i 2 .
(3.2)

In order to validate these asymptotic expressions we compare their predictions for the solitonic gas described in Section 3. The skewness ς and kurtosis κ are computed from numerical simulations and the time average is then taken. The results of the comparison are provided in ρ → 0. However, as ρ → 0, the jump from the initial value of ς(0) or κ(0) to the stationary one will diminish, since the interactions between solitons become sparser.

The asymptotic formulas (3.1), (3.2) can be used to estimate the skewness ς and kurtosis κ for various values of model parameters δ and S in initial stages of the solitonic gas evolution. However, these formulas contain the statistical averages. The proposed Monte-Carlo approach is briefly summarized here. Namely, we generate M random independent realizations of a solitonic gas, each sample consisting of N s solitons. The numerical values of parameters M and N s used in simulations are given in Table 1. Then we use the assumption of the rarefied gas along with the knowledge of the analytical soliton shape in order to estimate some statistical quantities of the gas (see equations (3.1), (3.2), for example). A big number of Monte-Carlo realisations allows to annihilate local fluctuations and to obtain robust estimations.

The numerical results are presented on Figures 9 -11 where the shadowed areas show the statistical error due to Monte-Carlo simulations (±σ, ±1.96σ, where σ is the estimated variance). From these results one can clearly see that for the fixed density ρ an increase in the Stokes-Ursell number S leads to an increase in ς and κ (see Figure 9). The BBM coefficient δ has the antagonistic effect as it is illustrated on Figure 10. When both parameters S and δ are increased simultaneously2 , the dependence of statistical quantities is not monotonic anymore (see Figure 11).

Conclusions and perspectives

In this study we presented several numerical experiments on the solitonic gas turbulence in the framework of an integrable KdV and a nonintegrable regularized KdV-BBM equation. The numerical results reported above generalize previous investigations [START_REF] Pelinovsky | Numerical modeling of the KdV random wave field[END_REF][START_REF] Dutykh | Geometric numerical schemes for the KdV equation[END_REF] where only a limited number (a few dozens) of solitons were used to represent a solitonic gas. Consequently, we reduce the statistical error according to the law of large numbers.

First of all, we showed that the probability distribution for the solitonic gas remains quasi-invariant during the system evolution for both KdV and KdV-BBM cases. The special attention was payed to the statistical characteristics such as kurtosis and skewness which measure the 'heaviness' of tails and the asymmetry of the free surface elevation distribution. In particular, using the asymptotic methods and Monte-Carlo simulations we showed that both skewness and kurtosis increase with the Stokes-Ursell number S and decrease when the BBM term coefficient δ. When both parameters S and δ are increased gradually and simultaneously, these effects are in competition: first we observe the increase of these statistical characteristics, but then, this tendency is inversed and they decrease after reaching their respective maximal values (see Figure 11). We would like to underline that the proposed Monte-Carlo methodology is much less computationally expensive than direct numerical simulations. Despite the small number of Monte-Carlo runs (M = 100) the estimated statistical error is sufficiently small for the purposes of this study. On the other hand, this approach is restricted, strictly speaking, to the situations where the solitons are well separated, since it does not take into account superpositions of solitons.

The present study opens a number of perspectives for future investigations. More general nonlinearities could be included into the model along with some weak dissipative and forcing effects. This could allow us to observe Kolmogorov spectra of a solitonic gas [START_REF] Zakharov | Kolmogorov Spectra of Turbulence I[END_REF]. The nonintegrable effects need some time to be accumulated. Consequently, even longer simulation times are needed. The interaction of a solitonic gas with a random radiation field has to be studied as well.

xFigure 1 .

 1 Figure 1. Solitary wave shape dependence on various parameters: (a) amplitude a = 1.0, 1.2, 1.5 for S = 1.0, δ = 0 (b) δ = 0.0, 1.0, 2.0 for a = 1.0, S = 1.0 (c) S = 1.0, 2.0, 3.0 for a = 1.0, S = 1.0 (d) simultaneous change of (δ, S) = (0.0, 1.0), (1.0, 2.0), (2.0, 3.0) for a = 1.0.

Figure 2 .

 2 Figure 2. Conservation of the invariant I 1 (t) during the KdV (a) and KdV-BBM (b) simulations.

Figure 3 .

 3 Figure3. The initial condition (a) and the final state (b) of a random solitonic gas simulated using the KdV equation (δ = 0). The final state of the KdV-BBM simulation is shown on panel (c). Please, note that the final simulation times T 1 = T 2 (see Table1for the values of T 1,2 ). Parameter τ 0 denotes the time needed for an average soltion to go over the whole computational domain.

Figure 4 .Figure 6 .

 46 Figure 4. Space-time plot of a random solitonic gas under the KdV (a) and KdV-BBM (b) dynamics. The time arrow is directed upwards (the initial state corresponds to the bottom line). The rectangular box in the upper image (a) shows the area zoomed in the next Figure 5.

Figure 7 .

 7 Figure 7. Probability distributions of the squared normalized free surface elevation.

Figure 9 .

 9 Figure 9. Initial kurtosis (left) and skewness (right) dependence on the Stokes-Ursell number S.

Figure 10 .

 10 Figure 10. Initial kurtosis (left) and skewness (right) dependence on the BBMterm coefficient δ.

Figure 11 .

 11 Figure 11. Initial kurtosis (left) and skewness (right) behaviour when one increases simultaneously δ and S.

  The solitary wave collisions in this equation were studied earlier by[START_REF] Francius | Wave dynamics in nonlinear media with two dispersionless limits for long and short waves[END_REF] [START_REF] Francius | Wave dynamics in nonlinear media with two dispersionless limits for long and short waves[END_REF] and[START_REF] Kalisch | Solitary wave collisions in the regularized long wave equation[END_REF] [21]

Table 2 .

 2 One can see that the simple analytical model presented in this Section is clearly able to describe rarefied solitonic gases. Moreover, one can infer from (3.1), (3.2) the behaviour of the skewness and kurtosis as the density of the gas decreases

	Skewness, ς	Kurtosis, κ
	KdV (δ = 0) 2.6111/2.6365 7.5402/7.7047
	KdV-BBM (δ = 2) 2.5482/2.5732 7.1759/7.3359

Table 2 .

 2 Comparison of the numerical (left) to the approximate analytical (right) estimations for the skewness (3.1) and kurtosis (3.2) on the data analysed in the previous section.

Sometimes it is also called the solitonic turbulence, e.g. in[START_REF] Zakharov | Soliton turbulence[END_REF][START_REF] Dyachenko | Soliton turbulence in nonintegrable wave systems[END_REF].

It means that we introduce an auxiliary homotopy parameter λ ∈ [0, 1] which parametrizes the simultaneous change of parameters S and δ in the following way:S = S min (1 -λ) + S max λ, δ = δ min (1 -λ) + δ max λ.In the computations presented below we took the values: S min = 0.1, S max = 6.0, δ min = 0.0 and δ max = 4.0.

Acknowledgments

Acknowledgments D. Dutykh acknowledges the support from ERC under the research project ERC-2011-AdG 290562-MULTIWAVE. E. Pelinovsky would like to thank for the support the Volk-swagenStiftung, RFBR grants (14-05-00092). The authors would like to thank Gennady El for stimulating discussions on the topics of integrability and solitons theory.

is too weak to modify significantly the probability distribution. Moreover, for comparison we plot also on the same Figure the standard normal (Gaussian) distribution. One can see that numerical results show much heavier tails than the standard distribution depicted on the same plot.

The probability distribution can be characterized by several parameters. Perhaps two most important characteristics are listed hereinbelow:

, which measures the heaviness of the spectrum tail.

, which measures the asymmetry of the spectrum with respect to the mean.

These quantities are defined in terms of the normalized free surface elevation moments µ n := η n 0 . We note that κ = 3 and ς = 0 for the normal (Gaussian) distribution N(0, 1). The evolution of these quantities is shown on Figure 8. The kurtosis κ is shown on top panels (a,b) and the skewness ς on the bottom (c,d). The KdV simulation results are represented on the left images (a,d) and the KdV-BBM on the right (b,d). One can see that the qualitative behaviour of these quantities is similar in integrable and nonintegrable cases. After a rapid initial transient period both quantities κ and ς enter in a quasistationary regime which consists of fast small amplitude (±1.6%) oscillations around the mean value. Initial values of the moments are higher, but then it drops down quickly due to soliton interactions. We note also that corresponding values of κ and ς are slightly lower in the KdV-BBM case due to the differences in solitary wave shapes.


Francis Filbet
email: filbet@math.univ-lyon1.fr

Lorenzo Pareschi
email: lorenzo.pareschi@unife.it

Thomas Rey
email: trey@cscamm.umd.edu

ON STEADY-STATE PRESERVING SPECTRAL METHODS FOR HOMOGENEOUS BOLTZMANN EQUATIONS

Keywords: 2010 Mathematics Subject Classification. Primary: 76P05, 82C40, Secondary: 65N35 Boltzmann equation, spectral method, Fourier series, steady-states preserving, micromacro decomposition, Maxwellian

In this note, we present a general way to construct spectral methods for the collision operator of the Boltzmann equation which preserves exactly the Maxwellian steady-state of the system. We show that the resulting method is able to approximate with spectral accuracy the solution uniformly in time.

Résumé. Dans cette note, nous présentons une construction générale de méthodes spectrales pour l'opérateur de collision de l'équation de Boltzmann permettant de préserver exactement les états stationnaires Maxwellien de ce type d'équations. Cette nouvelle approche est basée sur une décomposition de type "micro-macro" de la solution de l'équation, tout en restant très proche d'une méthode spectrale plus classique. Nous montrons que les méthodes obtenues sont capables d'approcher avec une précision spectrale, uniformément en temps, la solution de l'équation considérée, et nous présentons leur efficacité dans un test numérique.

1. Introduction Spectral methods for the numerical solution of the homogeneous Boltzmann equation have been proposed originally in [START_REF] Pareschi | A spectral method for the homogeneous Boltzmann equation[END_REF] and [START_REF] Pareschi | Numerical solution of the Boltzmann equation I. Spectrally accurate approximation of the collision operator[END_REF]. Recently their stability properties were studied in [START_REF] Filbet | Analysis of spectral methods for the homogeneous Boltzmann equation[END_REF]. Related approaches, based on the use of the Fourier transform have been introduced by [START_REF] Bobylev | Fast deterministic method of solving the Boltzmann equation for hard spheres[END_REF] and [START_REF] Gamba | Spectral-Lagrangian methods for collisional models of non-equilibrium statistical states[END_REF].

In classical kinetic theory the behavior of a rarefied gas of particles is described by the Boltzmann equation [START_REF] Cercignani | The Mathematical Theory of Dilute Gases[END_REF] ∂f ∂t

+ v • ∇ x f = Q(f, f ) where f (t, x, v), x ∈ R d , v ∈ R d (d ≥ 2)
, is the time-dependent particle distribution function in the phase space and the collision operator Q is defined by

(1) Q(f, f )(v) = R d S d-1 B(cos θ, |v -v * |) f ′ * f ′ -f * f dω dv * .
Time and position act only as parameters in Q and therefore will be omitted in its description.

In (1) we used the shorthands

f = f (v), f * = f (v * ), f ′ = f (v ′ ), f ′ * = f (v ′ *
). The velocities of the colliding pairs (v, v * ) and (v ′ , v ′ * ) are related by

v ′ = v + v * 2 + |v -v * | 2 σ, v ′ * = v + v * 2 - |v -v * | 2 σ.
The collision kernel B is a non-negative function which only depends on |v -

v * | and cos θ = ((v -v * )/|v -v * |) • ω.
In this note, we shall assume that the system considered is space homogeneous and then the particle distribution function f will only depends on the time t and the velocity v. The Boltzmann equation then becomes

(2) ∂f ∂t = Q(f, f ).
Boltzmann's collision operator has the fundamental properties of conserving mass, momentum and energy

v∈R d Q(f, f ) φ(v) dv = 0, φ(v) = 1, v, |v| 2
and satisfies the well-known Boltzmann's H-theorem

- d dt v∈R d f log f dv = - R d Q(f, f ) log(f ) dv ≥ 0.
The functional -f log f is the entropy of the solution. Boltzmann H-theorem implies that any equilibrium distribution function has the form of a locally Maxwellian distribution

M (ρ, u, T )(v) = ρ (2πT ) d/2 exp - |u -v| 2 2T ,
where ρ, u, T are the density, mean velocity and temperature of the gas:

ρ = R d f (v)dv, u = 1 ρ R d vf (v)dv, T = 1 dρ R d |u -v| 2 f (v)dv.
A major drawback of the classical spectral method is the lack of exact conservations and, as a consequence the incapacity of the scheme to preserve the Maxwellian steady states of the system. In this paper we overcome this drawback thanks to a new reformulation of the method which permits to preserve the spectral accuracy and to capture the long time behavior of the system.

A steady-state preserving spectral method

Although the method applies in principle to any collisional kinetic equation which possesses a local Maxwellian-type equilibrium, here we describe the method in the case of the Boltzmann operator.

A micro-macro decomposition. Let us start with the decomposition

(3) f = M + g,
with M the local Maxwellian equilibrium and g such that R d g φ dv = 0, φ = 1, v, |v| 2 . When inserted into a Boltzmann-type collision operator, the decomposition (3) gives ( 4)

Q(f, f ) = L(M, g) + Q(g, g),
where

L(M, g) = Q(g, M ) + Q(M, g
) is a linear operator and we used the fact that ( 5)

Q(M, M ) = 0.
There are two major features in the decomposition ( 4):

(1) it embeds the identity (5);

(2) the steady state of ( 4) is given by g = 0. To illustrate the method let us consider now the space homogenous equation (2) that we rewrite using the micro-macro decomposition as ( 6)

   ∂g ∂t = L(M, g) + Q(g, g), f = M + g.

Derivation of the spectral method.

To simplify notations we derived the spectral method in the classical setting introduced in [START_REF] Pareschi | Numerical solution of the Boltzmann equation I. Spectrally accurate approximation of the collision operator[END_REF], similarly it can be extended to the representation used in [START_REF] Mouhot | Fast methods for the Boltzmann collision integral[END_REF] for the derivation of fast algorithms. Thus, we perform the usual periodization in a bounded domain of the operators L and Q and denote by L λ and Q λ the operators with cut-off on the relative velocity on B 0 (2λπ).

Let us first set up the mathematical framework of our analysis. For any t ≥ 0, f N (v, t) is a trigonometric polynomial of degree N in v, i.e. f N (t) ∈ P N where

P N = span e ik•v | -N ≤ k j ≤ N, j = 1, . . . , d . Moreover, let P N : L 2 ([-π, π] 3 ) → P N be the orthogonal projection upon P N in the inner product of L 2 ([-π, π] 3 ) < f -P N f, φ >= 0, ∀ φ ∈ P N .
We denote the L 2 -norm by

||f || 2 = (< f, f >) 1/2 .
With this definition P N f = f N , where f N is the truncated Fourier series of f given by

f N (v) = N k=-N fk e ik•v , with fk = 1 (2π) d [-π,π] d f (v)e -ik•v dv.
We then write the Fourier-Galerkin approximation of the micro-macro equation (6)

   ∂g N ∂t = L λ N (M N , g N ) + Q λ N (g N , g N ), f N = M N + g N , where M N = P N M , g N = P N g, L λ N (M N , g N ) = P N L λ (M N , g N ) and Q λ N (g N , g N ) = P N Q λ (g N , g N ). (7)
More precisely, using multi-index notations in the above expressions we have

(8) Q λ N (g N , g N ) = N k=-N    N l,m=-N l+m=k ĝl ĝm β(l, m)    e ik•v , k = -N, . . . , N,
where the Boltzmann kernel modes β(l, m) = B(l, m) -B(m, m) are given by

(9) B(l, m) = B 0 (2λπ) S d-1
B(cos θ, |q|)e -i(l•q + +m•q -) dω dq.

In this last identity, q = v -v * is the relative velocity and q + , q -are the following parametrizations of the post-collisional velocities:

q + = 1 2 (q + |q|ω) , q -= 1 2 (q -|q|ω) .
It is immediate to show that Proposition 1. The function g N ≡ 0 is an admissible local equilibrium of the scheme [START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF] and therefore f N = M N is a local equilibrium state.

Spectral accuracy

In this section we show that ( 7) is a spectrally accurate approximation to [START_REF] Gamba | Spectral-Lagrangian methods for collisional models of non-equilibrium statistical states[END_REF] provided that M N is a spectrally accurate approximation of M . This is clearly guaranteed if initially the support is large enough. Note that due to space homogeneity, M does not change in time and so does M N . It is interesting to observe that the only difference between scheme (7) and the usual spectral method developed on the original formulation [START_REF] Mouhot | Fast methods for the Boltzmann collision integral[END_REF]

∂f N ∂t = ∂g N ∂t = Q λ N (f N , f N ), = L λ N (M N , g N ) + Q λ N (g N , g N ) + Q λ N (M N , M N ) is due to the constant (in time) term (11) Q λ N (M N , M N ) = 0
, which, as we will prove in the sequel, is spectrally small and is not present in [START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF].

If

f ∈ H r p ([-π, π] d ),
where r ≥ 0 is an integer and H r p ([-π, π] d ) is the subspace of the Sobolev space H r ([-π, π] d ), which consists of periodic functions, we have the following estimate for the spectral accuracy of the classical spectral method [START_REF] Pareschi | Numerical solution of the Boltzmann equation I. Spectrally accurate approximation of the collision operator[END_REF].

Theorem 1. Let f ∈ H r p ([-π, π] d ), r ≥ 0 then (12) ||Q λ (f, f ) -Q λ N (f N , f N )|| 2 ≤ C N r ||f || H r p + ||Q λ (f N , f N )|| H r p ,
Therefore we have the following result for ( 11)

Lemma 1. Let M N = PM where M is a local Maxwellian. Then (13) ||Q λ N (M N , M N )|| 2 ≤ C N r ||M || H r p + ||Q λ (M N , M N )|| H r p ,
Merging the two results we obtain the spectral accuracy of the new steady state formulation

Theorem 2. Let f ∈ H r p ([-π, π] 3 ), r ≥ 0 then ||Q λ (f, f ) -L λ N (M N , g N ) -Q λ N (g N , g N )|| 2 ≤ C N r ||f || H r p + ||M || H r p (14) +||Q λ (f N , f N )|| H r p + ||Q λ (M N , M N )|| H r p .

A numerical example

We present in this Section a numerical example of our method for the space homogeneous Boltzmann equation in dimension 2, with Maxwell molecules:

B(cos θ, |v -v * |) = 1 2π .
We compare the classical spectral method to the new steady-state preserving one. For this, we use an an exact solution of the homogeneous Boltzmann equation, the so called Bobylev-Krook-Wu solution [START_REF] Bobylev | Exact solutions of the Boltzmann equation[END_REF][START_REF] Krook | Exact solutions of the Boltzmann equation[END_REF]. It is given by

f BKW (t, v) = exp(-v 2 /2S) 2πS 2 2 S -1 + 1 -S 2 S v 2 with S = S(t) = 1 -exp(-t/8)/2.
For the resolution of the Boltzmann equation ( 2), we use the fast spectral method of [START_REF] Filbet | Solving the Boltzmann Equation in N log2 N[END_REF], with N = 32 half-modes in each direction of the box [-V, V ] 2 for V = 8. We take M = 8 angular discretizations. Here, we insists that both the classical and steady-states preserving methods rely on the same numerical algorithm. The only difference between the two is the presence of the constant in time term [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF]. In particular, both methods have the same computational cost, namely O (M N log 2 N ).

(f |M )(t) := R d f (t, v) log f (t, v) M (v) ,
and the absolute L 2 error between the numerical solution f (t, v) and the exact one f BKW (t, v). We observe that for both quantities, the behavior of the steady-state preserving method is better than the classical one. In particular, the new method achieves a nice monotonous decay of the relative entropy, without the large time increase of the classical spectral method. This is due to the fact that the equilibrium of this latter methods are constants [START_REF] Filbet | Analysis of spectral methods for the homogeneous Boltzmann equation[END_REF]. The behavior of the L 2 error is also monotonous for the steady-state preserving method, which is not the case for the classical one. We then present in Figure 2 a comparison between the temperatures of the solutions obtained with the classical and steady-state preserving methods. As expected with such a few Fourier modes, the temperature of the classical spectral method is not preserved, but grows almost linearly because of the convergence toward the discrete equilibrium, a constant on the cube. In contrast, it seems that the steady-state preserving spectral methods improve greatly this behavior: although not perfectly conserved, the temperature of the solution decays by roughly 0.1% and then stabilizes on a constant. Finally, as a last observation, we noticed that even if the positivity of the solution is not preserved by the new method, the number of nonpositive cells seems to be lower at a fixed time for the steady-state preserving method, compared to the classical one. We recall that it was already proved in [START_REF] Filbet | Analysis of spectral methods for the homogeneous Boltzmann equation[END_REF] that the number of nonpositive cells is asymptotically 0 for the classical spectral method: Theorem 3 (Filbet,Mouhot,[START_REF] Filbet | Analysis of spectral methods for the homogeneous Boltzmann equation[END_REF]). Let f N be the unique, global solution to the truncated, Fourier transformed Boltzmann equation, associated to a nonnegative initial datum f 0 ∈ H k p ([-π, π] d ) with k > d/2. Then the mass of the negative values of f N can be made uniformly in time small as the number of Fourier modes N → ∞.

Figure 1 .

 1 Figure 1. Comparison of the classical spectral method with the steady-state preserving one: Relative entropy H(f |M )(t) (left) and L 2 error f (t)-f BKW (t) 2 (right).

Figure 1

 1 Figure 1 presents a comparison between both methods for the relative entropy of the solution f with respect to the global Maxwellian M H(f |M )(t) :=

Figure 2 .

 2 Figure 2. Comparison of the classical spectral method with the steady-state preserving one: Evolution of the temperature.


Tuan-Hung Vu

Catherine Olsson

Ivan Laptev

Aude Oliva

Josef Sivic

Predicting Actions from Static Scenes

Keywords: Action prediction, scene recognition, functional properties

Human actions naturally co-occur with scenes. In this work we aim to discover action-scene correlation for a large number of scene categories and to use such correlation for action prediction. Towards this goal, we collect a new SUN Action dataset with manual annotations of typical human actions for 397 scenes. We next discover action-scene associations and demonstrate that scene categories can be well identified from their associated actions. Using discovered associations, we address a new task of predicting human actions for images of static scenes. We evaluate prediction of 23 and 38 action classes for images of indoor and outdoor scenes respectively and show promising results. We also propose a new application of geo-localized action prediction and demonstrate ability of our method to automatically answer queries such as "Where is a good place for a picnic?" or "Can I cycle along this path?".

Introduction

Our environments, such as living rooms, cafes and offices, vary in objects and geometry, but also in actions that we usually do in these places (e.g., we typically work in offices and cook or eat in kitchens). As illustrated in Figure 1, scene types are, indeed, often correlated with specific sets of typical actions. The goal of this work is to explore such correlation and to develop algorithms able to answer questions such as "What are typical actions for a given scene?", "Where is a good place to have a picnic?" or "Can I cycle along this path?". Automatic answers to such questions could be useful for several purposes. First, action prediction could provide scene-specific priors when recognizing human actions. For example, relaxing is common on beaches but not on streets; cooking is common in kitchens but not in offices. Second, deviations from an expected set of actions could be used to identify abnormal activities. Third, as we show in this paper, automatic action prediction for geo-localized images could support the search of places suited for particular purposes.

Computer vision has a rich body of work on recognizing human actions [START_REF] Laptev | Learning realistic human actions from movies[END_REF][START_REF] Niebles | Modeling temporal structure of decomposable motion segments for activity classification[END_REF][START_REF] Sadanand | Action bank: A high-level representation of activity in video[END_REF][START_REF] Wang | Action Recognition by Dense Trajectories[END_REF][START_REF] Marszalek | Actions in context[END_REF] and scenes [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF][START_REF] Oliva | Modeling the shape of the scene: A holistic representation of the spatial envelope[END_REF][START_REF] Quattoni | Recognizing indoor scenes[END_REF][START_REF] Vogel | Natural scene retrieval based on a semantic modeling step[END_REF]. Most of this work addresses the problems of action and scene recognition separately. Recently, several methods have shown advantages of recognizing actions or tracking people in the context of their environments [5, Fig. 1: Images of scene classes sandbar and temple east asia from the SUN dataset [START_REF] Jianxiong | Sun database: Largescale scene recognition from abbey to zoo[END_REF] together with probabilities for the five most likely actions, predicted manually by people (red) and by our method (green). 10]. Similarly, the interplay between human poses and objects has been studied in [START_REF] Grabner | What makes a chair a chair?[END_REF][START_REF] Gupta | From 3d scene geometry to human workspace[END_REF][START_REF] Delaitre | Scene semantics from long-term observation of people[END_REF]. While previous work has looked at functional properties for a few selected classes of scenes and objects, here we aim to exploit correlation between scenes and actions at a large scale of hundreds of scene categories. Using the discovered correlations, we demonstrate prediction of human actions for test images of outdoor scenes such as, for example, found on Google maps.

To reach our goal, we construct a new SUN Action dataset and collect manual annotations of human actions for 7940 images of 397 scene categories from the SUN dataset [START_REF] Jianxiong | Sun database: Largescale scene recognition from abbey to zoo[END_REF]. Analysis of this data reveals strong action-scene correlation for the majority of scene categories. Notably, we show that an image's scene category can be determined from corresponding textual descriptions of characteristic actions for that image.

Using the discovered relations between scenes and actions, we next address the task of automatic action prediction for images of static scenes. We consider 38 outdoor and 23 indoor action classes and associate typical action labels with 397 scene categories. Using such scene-based action annotation we learn visual classifiers for each action category and predict actions for images of static scenes as illustrated in Figure 1.

We finally demonstrate an application of our method to the new task of geo-localized action prediction. Our motivation comes from the large amount of publically-available geo-tagged images (e.g., on Flickr and Panoramio) which is expected to grow even faster with the introduction of new wearable devices such as Google Glass. Application of automatic action prediction to such images will enable the search for places based on their function, including specific actions such as swimming, having picnic, hiking and many others. In our experiments we use geo-localized images of outdoor scenes collected from panoramio.com and demonstrate examples of successful action prediction on the map of France.

In summary, we make the following three contributions. First, we present a new dataset with manual annotations of typical actions for 397 scene classes (see Section 3) and use it to analyze action-scene correlations (see Section 4). Second, based on the discovered correlations, we demonstrate successful action prediction for images of static scenes (see Section 5). Finally, we propose a new task of geolocalized action prediction. We apply our method to geo-tagged images on the web and show encouraging results of searching maps for locations suitable for particular activities (see Section 6).

Related work

Relatively few papers explore relations between scenes and actions. Li et al. [START_REF] Li | What, where and who? classifying events by scene and object recognition[END_REF] propose a graphical model combining evidence from object and scene categories for action recognition in still images. Marszalek et al. [START_REF] Marszalek | Actions in context[END_REF] propose a joint framework for scene and action recognition in video. While most of the work in action recognition targets actions depicted in images or video, here we address a different task and predict actions in scene images with no action observations. Action prediction has been recently addressed by Kitani et al. [START_REF] Kitani | Activity forecasting[END_REF] and Walker et al. [START_REF] Walker | Patch to the future: Unsupervised visual prediction[END_REF] aiming to model future motion of people and cars using priors derived from the scene. Yuen and Torralba [START_REF] Yuen | A data-driven approach for event prediction[END_REF] predict motion for images of static scenes by searching and transferring motion cues from video scenes with similar appearance. Our work complements these efforts and investigates action prediction for a large set of scenes and actions.

Recognition of functional properties of objects and scenes is an interesting but less explored area of computer vision. Relations between people and objects as well as between human poses and scene geometry have been investigated in [START_REF] Grabner | What makes a chair a chair?[END_REF][START_REF] Gupta | From 3d scene geometry to human workspace[END_REF][START_REF] Delaitre | Scene semantics from long-term observation of people[END_REF]. Patterson and Hays [START_REF] Patterson | Sun attribute database: Discovering, annotating, and recognizing scene attributes[END_REF] annotate scene images with a set of global attributes of various types (i.e: material, surface property, affordance and spatial envelope), and recognize attributes from scene images. Unlike any previous work, we here aim to model functional properties for a wide range of scene classes. Our work is similar in spirit to Arietta et al. [START_REF] Arietta | On relating visual elements to city statistics[END_REF] and Khosla et al. [START_REF] Khosla | Looking beyond the visible scene[END_REF] who aim to predict non-observed scene properties such as crime rate in the area.

Dataset

Dataset annotation. To analyze correlations over a wide range of scene categories and a rich set of actions, we gather the novel SUN Action dataset (short for "Scene UNderstanding -Action") with manual annotations of typical human actions for images of static scenes. We use scene images from the SUN dataset [START_REF] Jianxiong | Sun database: Largescale scene recognition from abbey to zoo[END_REF]. For each of the 397 well-sampled scene categories we collect freeform annotations of typical actions for the twenty "most typical" images in that category [START_REF] Ehinger | Estimating scene typicality from human ratings and image features[END_REF], for a total of 7940 images. Annotations were crowdsourced using Amazon Mechanical Turk (AMT) 1 . AMT workers were shown images of scenes and were asked to list between one and three words or short phrases for each scene describing a typical action that one would usually do there. Scene category labels were not provided. low response diversity). Such scenes often correspond to places that have been designed for specific purposes (tennis court) or where the natural environment limits the set of possible actions (wave). In contrast, scene classes with high entropy of responses (Figure 2, bottom) are places that afford many actions (e.g., a television studio, where many actions need to take place over the course of filming) or unfamiliar places (an anechoic chamber).

The majority of images in the SUN Action dataset contain no people. We found this property to be important for collecting unbiased annotations of typical actions. For a few images containing people we have observed a bias in action annotations towards actions depicted in the image. An example of such a bias is shown in Figure 3 illustrating two crosswalk scenes, one without people and one with a cycling person. In the scene containing the cyclist, the predominant response was "bike", unlike other images in the crosswalk category.

Processing of action responses. Action responses were gathered in free-form natural language and require preprocessing for our further analysis. Many of responses contain nearly identical information but differ in grammatical structure, such as "read the book while on the flight" and "read a book". Our first pass of preprocessing converts responses into simplified action annotations by extracting verbs or verb-noun patterns from each response. This strategy reduces the response space while preserving meaning. For example, responses like "read the book while on the flight" or "avoid eye contact with neighbours" are trimmed to "read book" and "avoid eye contact" respectively. We use the Stanford NLP toolbox [START_REF] Toutanova | Feature-rich part-of-speech tagging with a cyclic dependency network[END_REF] for part-of-speech tagging, stemming, and removal of stop words, and extract either verbs or verb-noun patterns from each response. Responses containing no verbs are removed. The words extracted in this stage of preprocessing are used as input to predict scene categories in Section 4.

For the action prediction task in Section 5 we manually group semantically similar action responses into action classes. To define action classes, we automatically extract 100 most frequent verb patterns, i.e. single verbs, verb+noun, etc., from action responses. Patterns with similar meaning are then manually merged yielding action labels, for example, "walk on grass" and "walk on sand" are merged into "walk". We note that the automatic parsing of natural language into action categories is an open problem beyond our work. In particular, we separate scenes into 197 outdoor and 203 indoor categories and define corresponding 38 outdoor and 23 indoor action classes as listed in Figure 7.

Given the average of 17.3 action responses per image in our database and a potentially large number of typical actions for a scene, our per-image annotation is not exhaustive. To address this problem, we assume that instances of the same scene category share the same functional properties. We found this assumption to be valid in most cases in our database. We therefore assign the same action labels to all instances of a given scene category using the following label propagation strategy. A scene category C is labeled by an action A if images of C are labeled with A at least 20 times. Following this procedure, for each action label A we obtain a set of positive scene categories. The negative scene categories for A are those containing no A labels for any of their images. Results of our preprocessing together with the original action responses are available from [23].

Analysis of scene-action correlation

Are different scene categories correlated with distinctive sets of actions? Scene categories are often defined by what you would typically do there: for example, in an office one would typically work, whereas in a kitchen one would typically cook. Indeed, most man-made scenes around us have been created to facilitate certain actions.

This section verifies and quantifies relations between actions and scenes. We demonstrate successful recognition of a large number of scene categories from associated actions descriptions. We further investigate the structure of actionscene correlations with a hierarchical clustering analysis.

Predicting scenes from actions

To verify the hypothesized correlation between scene categories and actions, we conduct two classification experiments using action annotations in the SUN Action dataset. We take inspiration from the field of text classification. In the SUN Action dataset, each image is associated with a collection of natural-language action descriptions. Classifying images based on a collection of associated responses is reminiscent of classifying documents based on their contents. However, there are two notable differences in our approach. First, the number of responses available per image (17.3 responses on average) is significantly lower compared to the number of words in a typical text document. Secondly, we wish to probe category membership using only a small collection of responses per image, to simulate asking a handful of people to provide a most typical action for the image and then performing classification based on the consensus of that set of responses. Therefore we use classification strategies that compare small queries to entire categories at a time.

Classification methods. We classify images using two simple bag-of-words techniques -Nearest Centroid and Naive Bayes. First, we divide the images in each class into 10 folds for cross-validation. Within the training set, the responses for each image are split into individual words. These word counts are combined and normalized across all images within a given class, to generate a word distribution histogram for each scene category. Within the test set, responses for each image are randomly grouped into chunks of 7 responses for that image, to simulate asking a handful of people at a time to provide a most typical action for each image. Responses within each chunk are then split into individual words to form bag-of-words queries.

In nearest centroid classification, the bag-of-words queries are normalized to form histograms, which are compared with category histograms according to histogram intersection distance. The scene category centroid with the smallest distance from the query is selected as the class label.

In Maximum Likelihood Naive Bayes classification, the category histograms are interpreted as empirical likelihood estimates: the likelihood Pr(w|c) of observing word w in association with an image of class c is assumed to be the The empirical likelihood estimate makes no explicit provision for estimating the likelihood of unobserved word-class pairs. To address this issue, we compute the minimum class-conditional likelihood over all words and classes in the dataset, min w,c (Pr(w|c)), and use this probability to stand in as the class-conditional likelihood for unobserved words. We assume a uniform prior over scene categories, enabling maximum likelihood estimation: that is to say, bag-of-words queries are classified according to which class provides the largest class-conditional likelihood.

Results. Figure 4 illustrates results of scene classification. As visualizing results across 397 individual classes is difficult, we select a 33-category subset of well-recognized and semantically important scene categories. To select the 33category subset, we have asked four of our collaborators to nominate 20-40 most important scene types. Out of the 80 scene types with most annotated images, 35 received at least two nominations and were slated for inclusion. "Cathedral" was removed for not being different enough from "church", and "abbey" and "coast" was removed for containing only aerial shots, leaving a final slate of 33 scene categories.

The confusion matrix for a Naive Bayes method in Figure 4(a) shows a strong diagonal indicating excellent classification performance. While most classes have almost perfect classification accuracy, a few classes are confused by the classifier due to the sharing of common actions. For example, scene categories "basement" and "attic" are both often annotated by actions "store" and "clean", while scene categories "river" and "lake" are frequently labeled with "swim" and "fishing".

Quantitative classification results of the two methods for the 33-category subset and all 397 scene classes are shown in Figure 4(b). Notably, both methods perform considerably better than chance while Naive Bayes provides better performance than Nearest Centroid. The fact that such simple classification methods yield very good performance indicates a strong correlation between scene categories and human actions: different scene categories have distinct patterns of associated actions. This confirms our initial hypothesis of a very strong relation between scene categories and their functional properties.

Action-based scene clustering

We seek to further investigate the structure of correlations between scene categories and actions: Which scene categories are more similar in terms of their function? We use hierarchical clustering and group scene descriptors at multiple scales. At the finest scale, only the most similar scene types cluster together, whereas at coarse scales, clusters are larger and encompass more dissimilar scene types. Dendrogram visualizations in Figure 5 show the progression of clustering patterns from fine to coarse: categories are represented as "leaves"; branchings closer to the leaves of the tree connect classes that cluster together under finegrained clusterings; and branchings closer to the trunk of the tree encompass broader clusterings. The height of each linkage in the dendrogram indicates the distance between the subclusters it connects.

The two dendrograms in Figure 5 illustrate image-based and action-based scene clustering. In the first case, distances between scenes were obtained as Eucledian distances of corresponding image descriptors (see Section 5). Clustering based on human action annotations was obtained using χ 2 -distances between scene representations in terms of bag-of-words histograms (see Section 4.1). We observe that image-based similarity in Figure 5(left) captures substantially different information about scene classes as compared to action similarity Figure 5(right). For example, "alley", "bathroom" and "subway interior" are grouped together according to visual similarity due to their similar geometrty and texture, but are separated according to action similarity since alleys, bathrooms and subways have different function. Another example is that visual similarity places "mountain, snowy" is a category of its own because no other class commonly depicts open white peaks, whereas action similarity places mountains together with other outdoor places that are associated with hiking, taking photos, and related actions.

Visual Action Prediction

People can easily determine appropriate actions to perform in a given place. Are machines able to do the same thing? In Section 4 we have addressed the related problem of predicting scene categories from a set of associated actions. Here we turn to the problem of predicting typical actions for an image of a scene. We approach the problem of visual action prediction using standard image classification techniques in terms of local features and binary classifiers. To train image classifiers we use action labels derived from action annotations as described in Section 3. We predict actions separately for indoor and outdoor scenes.

We test two different schemes for action prediction. Under the first scheme (S1), we train action classifiers directly from images using action labels only. Under the second scheme (S2), we first classify images into scene categories as an intermediate step and then assign action scores based on the obtained scores of scene classifiers. We assume that any particular test image belongs to one scene category only, therefore a score for an action a in a given image is defined as a max scene score over scene categories S associated with a.

Implementation Details

Image representation. Our image classification pipeline follows standard approaches and consists of densely extracted local image features, a learned visual vocabulary and a feature encoding step. For local image features we use HOG2x2 [START_REF] Jianxiong | Sun database: Largescale scene recognition from abbey to zoo[END_REF], SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] and CSIFT [START_REF] Abdel-Hakim | Csift: A sift descriptor with color invariant characteristics[END_REF][START_REF] Van De Sande | Empowering visual categorization with the gpu[END_REF] descriptors. Descriptor dimension is reduced by PCA. For the encoding phase we consider two popular encoding techniques: histogram encoding (BoW) and Fisher Vector encoding (FV). To exploit spatial information, we apply Spatial Pyramid framework [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF], using grids of size 1x1, 2x2, 3x1. Each grid cell is represented either by BoW or FV vectors. The resulting vectors are normalized and concatenated to create the final representation. In the rest of this section we use the format <descriptor > <encoding technique> to denote image representation techniques, e.g. CSIFT FV as Fisher Vector encoding for CSIFT descriptors.

Classification. For the classification, we train SVM classifiers using LIBSVM toolbox [START_REF] Chang | Libsvm: A library for support vector machines[END_REF]. Linear kernels are used for image representation by Fisher Vector. For the histogram representation (SIFT BoW /CSIFT BoW), we use χ 2 kernel [START_REF] Zhang | Local features and kernels for classification of texture and object categories: A comprehensive study[END_REF]. With HOG BoW, we exploit Histogram Intersection kernel. Training by SVM, we can boost up the performance by simply using a linear combination of kernels. In our experiments, we aggregate kernels with equal weights.

Experimental results

For training and testing the classifiers, we randomly divide the dataset into two equal parts. Our training and testing splits are balanced in number of images per scene category. Our results for action prediction are summarized in Figures 678and Table 1.

We use mean Average Precision (mAP) as the performance measure. To get mAP, we first compute the area under precision-recall curve, or Average Precision (AP), for each class. Then mAP is determined as the mean of average precisions across all classes. We obtain best prediction mAP of 60.99% for outdoor actions and 52.09% for indoor actions using combination of HOG BOW, HOG FV and CSIFT FV kernels. Our result is significantly higher than the mAP at chance level, i.e. 6.32% for outdoor action classes and 4.24% for indoor action classes. Figures 678are produced with our best kernel combination. In Figure 7, we show classification results with 38 outdoor and 23 indoor action classes sorted by AP. For better visualization of prediction results for some action classes, we show example images in Figure 8. The last two columns depict some hard positive and hard negative samples for each class. For outdoor scenes, action classes such as "hike", "pray", having rather typical color/structure, are easier to classify than other "can-do-almost-every-where" action classes like "learn". While people can often differentiate universities from other buildings based e.g., on the text and other cues, the task is still difficult for current vision systems, especially, for those exploiting global image representations. We notice that indoor actions are more structure-dependent than outdoor actions. In our experiment APs of indoor action classes are generally lower than APs of outdoor action classes. We also observe different levels of difficulty among indoor actions, e.g., detecting bowling lanes is easier than detecting sink-like structures. We found building action classifiers challenging, because positive samples are possibly images from very different scene categories, thus covering much larger range of visual texture and structure.

We also aggregate predicted action scores for test images and try to estimate the score contribution. Figure 6 shows some test images along with manual action annotations and automatic action predictions. For this visualization, we map SVM scores of test images to probabilities using Platt's sigmoid [START_REF] Platt | Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods[END_REF], with parameters estimated during the training phase. Even though the results are not perfect, we still observe a good match between distributions of annotated and predicted actions. Our predictors successfully give reasonable responses like "swim", "take a walk" and "relax" to a beach image, or "take picture", "pray" and "sightsee" to a temple image. Other qualitative results of action prediction by our method are available at [23].

For more quantitative analysis, we now consider Table 1. The table shows action prediction mAPs of two proposed training schemes combined with different image representation techniques. By comparing results of the two schemes, as shown in two columns S1 and S2, we can conclude that learning action classifiers directly achieve better prediction performance than aggregating mul- classification [START_REF] Perronnin | Fisher kernels on visual vocabularies for image categorization[END_REF]. Significant performance difference between SIFT and CSIFT proves that color information is useful for the task. Also, linear combination of multiple kernels does improve the performance. Among our three tested kernel combinations, using HOG BoW, HOG FV and CSIFT FV yields the best result.

In conclusion, we have shown high accuracy for a new task of action prediction evaluated on a large number of action and scene classes.

Image-based Geo-Mapping of Actions

One possible application of scene-based action prediction is to search for places in which to do a specific action. For example, a user may ask "Where can I camp in the Mont Blanc valley?" or "Where can I sunbathe in Tuscany?". Such queries are currently not supported by map services such as Google Maps or Bing Maps. To address this problem, we introduce Image-based Geo-Mapping of Action (IGMA), an application for geo-localizing actions on a map and answering map queries of the type "Where can I do X?". Results are derived from geo-localized scene images publicly available on the Internet. IGMA is the first attempt to automatically answer geo-localized action queries. Our strategy of predicting actions at a broad spatial scale using geo-localized images enables us to go beyond manual location-action labels: for example, one can "have a picnic" not only at a designated picnic area, but also in a grassy countryside field.

Collecting the Panoramio dataset. We use Panoramio image sharing service [START_REF]Panoramio service[END_REF] to collect a dataset with geo-localized images. Like the SUN Action dataset, the images in Panoramio contain few to no people. The Panoramio service provides a REST API for selecting images: given a range of longitude and latitude values, Panoramio returns a JSON file of image properties including im- age URL and geographical position. For our experiment, we collected Panoramio images of France, with longitude from -5 • to 8 • and latitude from 41 • to 51 • . In total, our dataset contains over 38, 000 distinct geo-tagged images.

Dense map of actions. Our goal is to construct dense maps visualizing places where people would likely perform certain actions. We construct these maps by applying the scene-based action classifiers computed in Section 5.2 to the collected Panoramio dataset using the following procedure. For a given action, we compute the top-scored Panoramio images. We generate a dense map from this list of scores and geo-locations by modeling the map using a Gaussian Mixture Model (GMM) with mixture components centered at the image locations and their weights set to corresponding action scores. The standard deviation σ for each component is set to a fixed value.

This initial dense map estimate is adjusted to compensate for non-uniform sampling of Panoramio images. Different population densities of the examined regions may introduce bias to the action density estimation. Therefore, we estimate the sampling density of Panoramio images. To estimate sampling density, we use the same GMM model above with the same σ for each Gaussian compo-nent as before, but with a uniform weight across all mixture model components rather than an action-score-based weight. The initial action map estimated from the highest scored images is normalized by the estimated sampling density of Panoramio images to correct the sampling density bias. We then get the final estimation of action density.

Figure 9 illustrates IGMA's suggestions for the question "Where can I ski in France?". We compare the estimated dense map produced by IGMA for the action "ski" with the the map of official skiing stations in France, acquired from [START_REF]Map of ski stations in france[END_REF]. Visually, our predictions have a high degree of correspondence with locations containing official skiing areas. Similarly, in Figure 10(left) we illustrate predictions for "hike" and "swim" in France. These results visually correspond to the sea-coast and mountain areas of France, confirming good geo-localization of actions.

Figure 10(right) illustrates an interesting result of predicting the "swim" action in Paris. This result suggests an area for further investigation: the recommended locations for swimming in Paris fall mainly along the river Seine, where swimming is very uncommon. While it is true that scene categories often have strong correlation with associated actions, not all scenes within a scene category share the same action affordances in practice. One possible approach to this issue might be to subdivide scene categories according to more fine-grained functional affordances, e.g. separately identifying rivers where people can and cannot swim.

Conclusion

In this work we have addressed a new problem of action prediction for a wide range of scene images. We have collected a new SUN Action dataset with manual annotations of typical actions for scene images, and discovered strong actionscene correlation for the majority of scene classes. Based on this correlation, we have learned to predict typical actions for a large set of scenes. Using standard state-of-the-art image classification techniques we have shown high accuracy of action prediction, which is an encouraging result for a new problem. To demonstrate potential advantages of our work, we have shown promising results on a new application Geo-Mapping of Actions (IGMA) enabling automatic answers to queries such as "Where can I do X?".

Fig. 2 :Fig. 3 :

 23 Fig. 2: SUN Action scene categories with corresponding histograms of action responses.Top row: Scene categories with low entropy of response histograms. Bottom row: scene categories with high entropy of response histograms. Low-entropy categories are often places designed for specific purposes (tennis court) or where the environment limits possible actions (wave). By comparison, high-entropy categories are places that afford many actions (television studio) or are unfamiliar (anechoic chamber)

Fig. 4 :

 4 Fig. 4: Results of action-based scene classification. (a): Confusion matrix for the 33category subset using Maximum Likelihood Naive Bayes estimation. The strong red line along the diagonal indicates excellent classification performance. A few pairs of categories e.g., (basement,attic) and (river,lake) are confused due to similarity in their characteristic actions. (b): Average accuracy (%) of scene classification for the 33category subset and for all 397 scene categories.

Fig. 5 :

 5 Fig.5: Results of hierarchical clustering of 33 scene categories based on the similarity of image descriptors (left) and action similarity (right). Image-based similarity groups similar-looking scenes despite their large difference in semantics such as "alley" and "bathroom". In contrast, action-based similarity results in more semantically meaningful clusters. For example, "mountain, snowy" is placed in a category of its own according to the visual similarity, whereas it is grouped together with other outdoor places on the basis of action similarity. number or observations of w within the class c responses in the training set divided by the total number of words in all class c responses combined. The word observation likelihoods are assumed to be conditionally independent (the "naive Bayes assumption"), enabling us to compute the class-conditional likelihood of each bag-of-word query as the product of each constituent word's empirical class likelihood: Pr(w 1 , w 2 , . . . w n |c) = Pr(w 1 |c) × Pr(w 2 |c) × . . . × Pr(w n |c).The empirical likelihood estimate makes no explicit provision for estimating the likelihood of unobserved word-class pairs. To address this issue, we compute the minimum class-conditional likelihood over all words and classes in the dataset, min w,c (Pr(w|c)), and use this probability to stand in as the class-conditional likelihood for unobserved words. We assume a uniform prior over scene categories, enabling maximum likelihood estimation: that is to say, bag-of-words queries are classified according to which class provides the largest class-conditional likelihood.

Fig. 6 :

 6 Fig. 6: Automatic visual action prediction for test images in SUN Action dataset.

Fig. 7 :Fig. 8 :

 78 Fig. 7: Results of action prediction for all 38 outdoor actions (top) and 23 indoor actions (bottom) sorted in the decreasing order of Average Precision (AP).Action Precision Recall

Fig. 9 :Fig. 10 :

 910 Fig. 9: "Where can I ski in France?" -(Top left) Official skiing stations in France [32]. (Top middle) Suggested places for skiing by IGMA. (Top right) Dense map of action "ski" generated by IGMA. (Bottom) Panoramio images of suggested places for skiing. hike -France swim -France swim -Paris

Table 1 :

 1 Scene-based outdoor (cyan) and indor (orange) action prediction results with different approaches. Note that mAP at chance level is 6.32% (outdoor) an d 4.24% (indoor). S1 and S2 columns respectively show classification mAP (%) of the two aforementioned training schemes. Column (S1-S2) shows the different mAP between two schemes. We observe consistently better performance of scheme S1: directly training binary action classifiers over scheme S2: aggregating scene classifiers.

	Method	S1	S2 S1-S2 S1	S2 S1-S2
	SIFT BoW	40.92 40.68 0.24 31.75 28.71 3.04
	SIFT FV	41.15 34.51 6.64 31.04 27.13 3.91
	CSIFT BoW	47.78 44.43 3.35 32.53 27.90 4.63
	CSIFT FV	49.52 41.65 7.87 36.29 29.70 6.59
	HOG BoW	47.03 45.93 1.10 37.91 35.50 2.41
	HOG FV	52.66 47.75 4.91 42.78 43.89 -1.11
	HOG BoW+ SIFT FV+ CSIFT FV 56.60 50.06 6.54 46.11 40.04 6.07
	HOG BoW+ HOG FV+ CSIFT FV 60.99 54.25 6.74 52.09 45.98 6.11
	SIFT FV+ HOG FV+ CSIFT FV 56.48 49.61 6.87 45.76 41.41 4.35

AMT workers gave consent (set by the MIT IRB) for each HIT they chose to perform.

Acknowledgements

This work is partly funded by ERC Activia, US National Science Foundation grant 1016862, Google Research Awards, MSR-INRIA laboratory and EIT-ICT labs.