Thierry E Huillet
email: thierry.huillet@cyu.fr

On random population growth punctuated by geometric catastrophic events

Keywords: Markov chain, random population growth; geometric catastrophe events, recurrence/transience transition, height and length of excursions, extinction events, time to extinction, large deviations PACS 02.50.Ey, 87.23. -n MSC primary 60J10, secondary 92A15

Catastrophe Markov chain population models have received a lot of attention in the recent past. Besides systematic random immigration events promoting growth, we study a particular case of populations simultaneously subject to the effect of geometric catastrophes that cause recurrent mass removal. We describe the subtle balance between the two such contradictory effects.

Introduction

Besides systematic immigration events, populations are often subject to the effect of catastrophic events that cause recurrent mass removal. This results in a subtle balance between the two contradictory effects. Many Markov chains models have been designed in an attempt to explain the transient and large-time behavior of the size of such populations. Some results concern the evaluation of the risk of extinction and the distribution of the population size in the case of total disasters where all individuals in the population are removed simultaneously, ( [START_REF] Swift | Transient probabilities for a simple birth-death-immigration process under the influence of total catastrophes[END_REF], [START_REF] Goncalves | Scaling features of two special Markov chains involving total disasters[END_REF]). Other works deal with the basic immigration process subject either to binomial or geometric catastrophes; that is when the population size is iteratively reduced according either to a binomial or a geometric law ( [START_REF] Brockwell | Birth, immigration and catastrophe processes[END_REF], [START_REF] Neuts | An interesting random walk on the non-negative integers[END_REF]). Such Markov chains are random walks on the non-negative integers (as a semigroup) which differ from standard random walks on the integers (as a group) in that a one-step move down from some positive integer cannot take the walker to a negative state, resulting in transition probabilities being state-dependent.

The binomial effect is appropriate when, on a catastrophic event, the individuals of the current population each die or survive in an independent and even way, resulting in a drastic depletion of individuals at each step, ( [START_REF] Neuts | An interesting random walk on the non-negative integers[END_REF], [START_REF] Ben-Ari | A random walk with catastrophes[END_REF], [START_REF] Fontes | Metastability of a random walk with catastrophes[END_REF]).

The geometric effect corresponds to softer deletion cases where the sequential decline in the population keeps going on but is stopped as soon as the current pop-ulation size is exhausted. In more precise words, if a geometric catastrophic event occurs, given the chain is in state x, the population shrinks by a random (geometrically distributed) amount so long as this amount does not exceed x; if it does the population size is set to 0 (a disaster event). In [START_REF] Neuts | An interesting random walk on the non-negative integers[END_REF], this model was called for this reason the 'truncated geometric' catastrophe model. This may be appropriate for some forms of catastrophic epidemics or when the catastrophe has a sequential propagation effect like in predator-prey models: the predator kills preys until it becomes satisfied but so long as resources are available. More examples and motivations can be found in Artalejo et al. [START_REF] Artalejo | Evaluating growth measures in an immigration process subject to binomial and geometric catastrophes[END_REF] and in Cairns and Pollett [START_REF] Cairns | Extinction times for a general birth, death and catastrophe process[END_REF], essentially in continuous-times.

In this work, in the spirit of [START_REF] Neuts | An interesting random walk on the non-negative integers[END_REF], we revisit the latter geometric catastrophe model in discrete-time, as a Markov chain on the non-negative integers.

Using a generating function approach, we first discuss the condition under which this process is recurrent (either positive or null) or transient. The recurrence/transience transition is sharp.

-In the positive-recurrent case, we describe the shape and features of the invariant probability measure. We emphasize that in the null-recurrent case, no non trivial invariant measure exists. Focus is then put on the time spent by the walker in the ground state {0} . We then show how to compute the double (space/time) generating functional of the process. Using this representation, we derive the generating functions of the first return time to zero (the length of the excursions) and of the first extinction time when the process is started at x 0 > 0. In the recurrent case, extinction occurs with probability 1. The harmonic (or scale) function of a version of the process forced to be absorbed in state {0} is used to give an expression of the law of the height of the excursions (the maximum population size reached between any two consecutive extinctions).

-In the transient case, after a finite number of visits to {0}, the chain drifts to ∞. We first emphasize that no non trivial invariant measure exists either. Using the expression of the double generating functional of the process, we get access to a precise large deviation result. The harmonic (or scale) function of a version of the process forced to be absorbed in state {0} is then used to give an expression of the probability of a first extinction. The law of the height of the excursions is also derived, with an atom at ∞. More complex growth/collapse or decay/surge models in the same vein, although in the continuum, were considered in [START_REF] Eliazar | Growth-collapse and decay-surge evolutions, and geometric Langevin equations[END_REF], [START_REF] Eliazar | Nonlinear Shot Noise: From aggregate dynamics to maximal dynamics[END_REF] and [START_REF] Eliazar | The maximal process of nonlinear shot noise[END_REF], where physical applications are developed, such as queueing processes arising in the Physics of dams or stress release issues arising in the Physics of earthquakes.

2 Geometric catastrophes, [START_REF] Brockwell | Birth, immigration and catastrophe processes[END_REF], [START_REF] Neuts | An interesting random walk on the non-negative integers[END_REF] Let us introduce the process under concern, as a Markov chain.

The model. Birth (growth): Let (β n ) n≥1 be an independent identically distributed (iid) sequence taking values in N := {1, 2, ...} , with b x := P (β 1 = x), x ≥ 1. We shall let B (z) := E z β be the common probability generating function of the βs.

Death (depletion): Let (δ n ) n≥1 be an iid shifted geometric(α) -distributed sequence 1 , with success parameter α ∈ (0, 1). With α := 1 -α, δ 1 then takes values in

N 0 := {0, 1, 2, ...} with P (δ 1 = x) = d x = αα x , x ≥ 0.
Consider the Markov chain with temporal evolution, [START_REF] Neuts | An interesting random walk on the non-negative integers[END_REF]:

X n+1 = X n + β n+1 with probability p = (X n -δ n+1 ) + with probability q = 1 -p. (1)
At each step n, the walker moves up with probability p and the amplitude of the upward move is β n+1 . The number of step-wise removed individuals (whenever available) is δ n+1 , with probability q : the distribution of the sizes of the catastrophes δ n+1 does not depend on the current size of the population, so long as there is enough prey. Stated differently, given the population size is x at n, the magnitude of the downward jump is x ∧ δ n+1 . Note that if X n = 0, then X n+1 = β n+1 with probability p (reflection at 0) and X n+1 = 0 with probability q (absorption at 0).

The one-step stochastic transition matrix P (obeying P 1 = 1, where 1 is a column vector of ones) of the Markov chain {X n } is: P (0, 0) = q, P (0, y) = pb y , y ≥ 1 P (x, y) = qd x-y = qαα x-y , x ≥ 1 and 1 ≤ y ≤ x P (x, 0) = q k≥x d k = qα x , x ≥ 1 and 0 = y < x P (x, y) = pb y-x , x ≥ 1 and y > x.

(

With π n := (P x0 (X n = 0) , P x0 (X n = 1) , ...) the column vector2 of the states' occupation probabilities at time n, π n+1 = π n P , X 0 d ∼ δ x0 , is the master equation of its temporal evolution. An equivalent recurrence for the probability generating function of {X n } started at X 0 = x 0 is:

Lemma 1 With Φ n (z) := E x0 z Xn = x≥0 z x P x0 (X n = x) , such that Φ 0 (z) = z x0 , it holds, Φ n+1 (z) = pB (z) - qαz α -z Φ n (z) + q 1 + αz α -z Φ n (α) . (3)
Proof:

From (2), Φ n+1 (z) = π n+1 (0) + pB (z) Φ n (z) + q y≥1 z y x≥y π n (x) αα x-y = π n+1 (0) + pB (z) Φ n (z) + qαz α -z (Φ n (α) -Φ n (z)) Φ n+1 (0) = π n+1 (0) = qΦ n (α) . 2
Remark: The geometric Markov chain is time-homogeneous, irreducible and aperiodic. As a result, all states are either recurrent or transient.

The recurrent case

The geometric Markov chain being either transient or recurrent, we shall first exhibit when this transition occurs. We shall then deal first with the recurrent case before switching to the transient regime.

Existence and shape of an invariant probability measure

Theorem 2 [START_REF] Neuts | An interesting random walk on the non-negative integers[END_REF]. The chain {X n } with geometric catastrophe is ergodic if and only if ρ := E (β 1 ) < ρ c := qp -1 α α < ∞. The probability generating function of X ∞ is then given by (4).

Proof: If a limiting random variable X ∞ exists, from [START_REF] Brockwell | Birth, immigration and catastrophe processes[END_REF], it has probability generating function Φ ∞ (z

) := E z X∞ obeying Φ ∞ (z) = α (1 -z) qΦ ∞ (α) (α -z) (1 -pB (z)) + qαz =: N (z) D (z) . Suppose ρ := B (1) = E (β 1 ) < ∞.
Only in such a case does the numerator N and denominator D both tend to 0 as z → 1 while the ratio Φ ∞ tends to 1, as required for Φ ∞ (z) to be a probability generating function. By l'Hospital rule

N (z) D (z) → 1 as z → 1, leading to Φ ∞ (α) = qα -pαρ qα . Thus π (0) := P (X ∞ = 0) = Φ ∞ (0) = qΦ ∞ (α) = q -pαρ/α,
and, with Φ Y (z) := q/ (1 -pB (z)) the probability generating function of a random variable Y obtained as a shifted-geometric(q) convolution of the βs

Φ ∞ (z) = (1 -z) (qα -pαρ) (α -z) (1 -pB (z)) + qαz = Φ Y (z) (1 -z) (qα -pαρ) (α -z) q + qαzΦ Y (z) π (x) : = π ∞ (x) = [z x ] Φ ∞ (z) = Φ ∞ (0) [z x ] 1 -z (α -z) (1 -pB (z)) + qαz , x ≥ 1.

If and only if

ρ = E (β 1 ) < ∞ and ζ := α 1 + q -1 pρ < 1 (ρ < ρ c ), then π (0) = q α (1 -ζ) ∈ (0, 1)
and the Markov chain is ergodic (positive recurrent and aperiodic).

The probability generating function of the invariant probability measure then takes the alternative form

Φ ∞ (z) = (1 -ζ) Φ Y (z) 1 -α 1 -zΦ Y (z) 1 -z -1 . Note that, with Φ Z (z) := 1-zΦ Y (z) (1+m)(1-z) = 1-Φ Y +1 (z) E(Y +1)(1-z)
the probability generating function of some random variable Z, this is also

Φ ∞ (z) = Φ Y (z) 1 -ζ 1 -ζΦ Z (z) , (4)
which is the product of the probability generating function of Y times the one of a shifted compound geometric(1 -ζ) random variable with compounding probability generating function Φ Z (z) (the delay probability generating function of

Y + 1). Note ζ < 1 ⇔ pE (β 1 ) = pρ < q α α = qE (δ 1 ) , observing α/α = E (δ 1 )
. This is a sub-criticality condition stating that the one-step unlimited average move down qE (δ 1 ) must exceed the average move up pρ of the geometric chain. 2

Corollary 3 If ρ < ρ c , the random variable X ∞ exists and is infinitely divisible (compound Poisson).

Proof: Consider first the random variable Y with probability generating function

Φ Y (z) := E z Y = q 1-pB(z) . It is a compound shifted-geometric random variable with compounding random variable β. It is infinitely divisible because Φ Y (z) = exp -λ (1 -ψ (z))
for some λ > 0 and some probability generating function ψ (z) with ψ (0) = 0. Indeed, with q =: e -λ , there exists a probability generating function

ψ solving Φ Y (z) := q/ (1 -pB (z)) = exp -λ (1 -ψ (z)) . It is ψ (z) = -log (1 -pB (z)) -log q ,
which is recognized as the probability generating function of a Fisher-log-series random variable. Shifted-geometric convolution random variables are infinitely divisible. For the same reason, the random variable with probability generating function

1-ζ 1-ζΦ Z (z)
is infinitely divisible as a shifted-compound geometric random variable with compounding random variable Z. From (4), Φ ∞ (z) is the product of these two probability generating functions and so is itself the probability generating function of an infinitely divisible random variable X ∞ . 2

Miscellaneous:

-Neuts in [START_REF] Neuts | An interesting random walk on the non-negative integers[END_REF] observes that if β is itself geometrically distributed (more generally of phase-type), then so are Y, Z, together with X ∞ , so with an explicit expression of π (x). The random variable X ∞ in particular admits geometric tails. Discrete phasetype random variables are the ones obtained as first hitting times of the absorbing state of a terminating Markov chain with finitely many states, [START_REF] Neuts | Matrix -Geometric-Solutions in Stochastic Models: an algorithmic approach[END_REF].

-In the general case for the law of β, an expression of π (x) := [z x ] Φ ∞ (z) can be obtained from (4), in principle, using the Faa di Bruno formula ( [START_REF] Comtet | Analyse combinatoire[END_REF], p. 146), observing Φ ∞ (z) = H (B (z)) for some generating function H, as a composition of generating functions.

-We also have: µ

:= E (X ∞ ) < ∞ ⇔ E Y 2 < ∞ ⇔ E β 2 < ∞, owing to: E Y 2 = 2E (Y ) 2 + E (G) E β 2 , E (G) = p/q. Specifically, if E β 2 < ∞ µ = m + α E Y 2 -E (Y ) 2 2 (1 -ζ) = m 1 + α E (β) E (β) m + E β 2 2 (1 -ζ) .
-More generally, it can also be checked that if β only has moments of order q < 1 + a (a > 0), then X ∞ only has moments of order q < a.

-

Finally, with m = Φ Y (1) = E (Y ) = q -1 pρ = q -1 pE (β), as z → 1 1 -zΦ Y (z) 1 -z ∼ 1 -z (1 -m (1 -z)) 1 -z ∼ 1 + m, so that Φ ∞ (1) ∼ (1 -ζ) [1 -α (1 + m)] -1 = 1.
-In the positive recurrent case, the geometric process is not time-reversible with respect to the invariant probability measure π: there is no detailed balance.

-Considerable simplifications are expected when β d ∼ δ 1 : in this case, the transition matrix P in (2) is of Hessenberg type, that is lower triangular with a non-zero upper diagonal. The corresponding process is skip-free to the right.

Time spent in the ground state {0}

Whenever the process {X n } is ergodic, it visits infinitely often all the states, in particular the state {0}, and a sample path of it is made of iid successive non-negative excursions (arches) through that state. By the ergodic theorem, the asymptotic fraction of time spent by {X n } in this state is

1 N N n=1 1 (X n = 0) → π (0) = P (X ∞ = 0) , almost surely, as N → ∞. (5)
Let τ 0,0 := inf (n ≥ 1 : X n = 0 | X 0 = 0) be the first return time to state {0} . By Kac's theorem, [START_REF] Kac | Random walk and the theory of Brownian motion[END_REF], its expected value is E (τ 0,0 ) = 1/π (0) where π (0

) = α 1 + q -1 pρ = q α (1 -ζ).
Suppose {X n } enters state {0} from above at some time n 1 . The first return time to state {0}, τ 0,0 := inf (n > n 1 : X n = 0 | X n1 = 0), is:

-either 1 if X n1 stays there with probability P (0, 0) = q in the next step; this corresponds to a trivial excursion of length 1 and height 0.

-or, with probability p = 1-q, {X n } starts a 'true' excursion with positive height and length τ + 0,0 ≥ 2. Thus,

E (τ 0,0 ) = 1 π (0) = 1 + pE τ + 0,0 and E τ + 0,0 = 1 p 1 π (0) -1 ≥ 2, (6)
entailing π (0) < 1/ (1 + 2p). Given {X n } enters state {0} from above at some time n 1 , it stays there with probability P (0, 0) = q in the next step, so {X n } will quit state {0} at time n 1 +G where G is a geometric random time with success probability p = 1 -q. After n 1 + G, the chain moves up before returning to state {0} again and the time it takes is τ + 0,0 . Considering two consecutive instants where {X n } enters state {0} from above (defining an alternating renewal process), the mean fraction of time θ spent in state {0} is:

θ := E (G) E (G) + E τ + 0,0 .
From the expression E (G) = 1/p and the value of E τ + 0,0 , consistently with (5), we get:

Proposition 4 In the positive recurrent case, the mean fraction of time θ spent in state {0} is:

θ = 1 1 + pE τ + 0,0 = π (0) . ( 7
)
3.3 No non trivial ( = 0) invariant measure in the null-recurrent case

If ζ = 1 (or ρ = ρ c ), the critical chain is null-recurrent with π (0) = qα -pαρ α = q α (1 -ζ) = 0 ⇒ π (x) = 0 for all x ≥ 1.
The chain has no non trivial ( = 0) invariant positive measure. It is not Harris-recurrent, [START_REF] Harris | The existence of stationary measures for certain Markov processes[END_REF], [START_REF] Harris | Transient Markov chains with stationary measures[END_REF] and [START_REF] Jamison | Markov chains recurrent in the sense of Harris[END_REF].

The generating functional of the geometric model

With x 0 ≥ 1, defining the double generating function

Φ x0 (u, z) = n≥0 u n Φ n (z) = n≥0 u n E x0 z Xn , from (3)
, we get

1 u (Φ x0 (u, z) -z x0 ) = pB (z) - qαz α -z Φ x0 (u, z) + q 1 + αz α -z Φ x0 (u, α) , together with Φ x0 (u, 0) = quΦ x0 (u, α) .
We obtain

Φ x0 (u, z) = z x0 (α -z) + qαuΦ x0 (u, α) (1 -z) (α -z) [1 -puB (z)] + qαzu (8) = z x0 (α -z) + α (1 -z) Φ x0 (u, 0) (α -z) [1 -puB (z)] + qαzu . So far, Φ x0 (u, z) is unknown since it requires the knowledge of Φ x0 (u, α) or Φ x0 (u, 0). Letting Φ x0 (u, z) =: N (u, z) /D (u, z), the denominator D (u, z) vanishes at u = u (z) = z -α qαz + p (z -α) B (z) . (9)
Note u (1) = 1 and u

(1) = 1 α (qα -αpρ) = q α (1 -ζ) . (10)
The generating function Φ x0 (u, z) is well-defined when u < u (z) and possibly when u = u (z) as in the recurrent case.

In the recurrent case (ζ ≤ 1), u (z) is concave and monotone increasing on the interval [α, 1] , owing to u (1) ≥ 0. The function u (z) has a well-defined inverse z (u) which maps [0, 1] to [α, 1]; this inverse is monotone increasing and convex on this interval. Because in the recurrent case, [u n ] Φ x0 (u, z) → Φ ∞ (z) as n → ∞, Φ x0 (u, z) also converges as z → z (u) . So both the numerator N and the denominator D of Φ x0 (u, z) must tend to 0 as z → z (u), meaning (by L'Hospital rule) that lim z→z(u)

N (u, z) D (u, z) = lim z→z(u) N (u, z) D (u, z) . Near z = z (u) , N (u, z) = A (z) + B (z) Φ x0 (u, 0) ∼ A (z (u)) + B (z (u)) Φ x0 (u, 0) + (z -z (u)) [A (z (u)) + B (z (u)) Φ x0 (u, 0)] D (u, z) = C (z) -D (z) u ∼ (z -z (u)) [C (z (u)) -D (z (u)) u] , imposing A (z (u)) + B (z (u)) Φ x0 (u, 0) = 0, thereby fixing Φ x0 (u, 0) = z (u) x0 (z (u) -α) α (1 -z (u)) =: G x0,0 (u) ≥ 0. ( 11
)
Note G x0,0 (1) = ∞ and G x0,0 (0) = P x0 (X 0 = 0) = 0. The function

G x0,0 (u) = n≥1 u n P x0 (X n = 0) = n≥1 u n P n (x 0 , 0) ,
is the Green kernel of the chain at the endpoints (x 0 , 0), [START_REF] Neveu | Chaines de Markov et théorie du potentiel[END_REF]. The matrix element P n (x 0 , 0) is the contact probability at 0 at time n, starting from x 0 . From ( 8) and ( 11), we thus get a closed form expression of Φ x0 (u, z) when x 0 ≥ 1, as (u ∈ [0, 1]) Proposition 5 In the recurrent case (ζ ≤ 1), with z (u) the inverse of u (z) defined in [START_REF] Eliazar | Growth-collapse and decay-surge evolutions, and geometric Langevin equations[END_REF],

Φ x0 (u, z) = z (u) x0 (z (u) -α) (1 -z) -z x0 (z -α) (1 -z (u)) (1 -z (u)) [(α -z) (1 -puB (z)) + qαzu] . (12)
Remark: If z = 1 Φ x0 (u, 1) = 1/ (1 -u) and Φ x0 (0, z) = z x0 , as required.

First return time to 0

When x 0 = 0, observing now G 0,0 (0

) = P 0 (X 0 = 0) = 1, G 0,0 (u) = 1+Φ 0 (u, 0) = 1+ z (u) -α α (1 -z (u)) = αz (u) α (1 -z (u))
, with G 0,0 (1) = ∞ (z (1) = 1).

(13) This function is the Green kernel at the endpoints (0, 0) .

If n ≥ 1, from the recurrence P 0 (X n = 0) =: P n (0, 0) = n m=0 P (τ 0,0 = m) P n-m (0, 0) , we see, from taking the generating function of both sides and observing the right hand-side is a convolution, that the pgf φ 0,0 (u) = E (u τ0,0 ) of the first return time to 0, τ 0,0 and G 0,0 (u) are related by the Feller relation (see [START_REF] Feller | An introduction to probability theory and its applications[END_REF] and [START_REF] Bingham | Stochastic Processes: Theory and Methods[END_REF] pp 3 -4 for example): G 0,0 (u) = 1+ G 0,0 (u) φ 0,0 (u). Hence, with φ 0,0 (0) = 0,

φ 0,0 (u) = E (u τ0,0 ) = 1 - 1 G 0,0 (u) = z (u) -α αz (u) , φ 0,0 (1) = 1. ( 14
)
In particular, observing z (1) = 1/u (1) = α q(1-ζ) , we get

E (τ 0,0 ) = αz (1) /α = α q (1 -ζ) if ζ < 1 (positive recurrence), = ∞ if ζ = 1 (null recurrence).
Note as required from Kac's theorem: E (τ 0,0 ) = 1/π (0) .

Contact probability at 0 and first time to extinction

Also, with x 0 ≥ 1 (φ x0,0 (0) = 0), using [START_REF] Eliazar | The maximal process of nonlinear shot noise[END_REF] and ( 13),

E (u τx 0 ,0 ) = φ x0,0 (u) = G x0,0 (u) G 0,0 (u) = 1 α z (u) x0-1 (z (u) -α) = φ 0,0 (u) z (u) x0 , ( 15
)
gives the probability generating function of the first hitting time of 0, starting from x 0 ≥ 1 (the first extinction time of the chain). We also get

E (τ x0,0 ) = z (1) α (αx 0 + α) = αx 0 + α q (1 -ζ) if ζ < 1 (positive recurrence), = ∞ if ζ = 1 (null recurrence). If z (1) < ∞ (σ 2 (β) < ∞)
, the variance of τ x0,0 in the positive-recurrent case is finite and found to be

σ 2 (τ x0,0 ) = z (1) + z (1) α (αx 0 + α) - z (1) 2 α 2 α 2 x 0 + 1 -α 2 = x 0 z (1) + z (1) -z (1) 2 + α α 2 α z (1) + z (1) -(1 + α) z (1) 2 . If σ 2 (β) = ∞, σ 2 (τ x0,0 ) = ∞.
We observe from ( 14)

φ 0,0 (u) = E (u τ0,0 ) = z (u) -α αz (u) , entailing z (u) = α 1 -αφ 0,0 (u)
.

The function φ 0,0 (u) being a probability generating function, z (u) is the probability generating function of some random variable τ ≥ 0; namely, the one of a compound shifted-geometric(α) probability generating function with compounding probability generating function φ 0,0 (u). And z (u) x0 is the probability generating function of the random variable

τ x0 = x0 k=1 τ k , ( 16
)
where τ k is an iid sequence with τ k ∼ τ . We have thus proven Theorem 6 In the recurrent case, the function z (u) is the probability generating function of some random variable τ ≥ 0. With τ x0 defined in [START_REF] Harris | Transient Markov chains with stationary measures[END_REF], the first time to extinction random variable τ x0,0 is decomposable as

τ x0,0 d = τ x0 + τ 0,0 , (17)
where τ x0 ≥ 0 and τ 0,0 ≥ 1 are independent. With φ 0,0 (u) given by ( 14), the probability generating function of τ x0,0 is given in [START_REF] Harris | The existence of stationary measures for certain Markov processes[END_REF] with z (u) = α/ (1 -αφ 0,0 (u)).

The function u (z) is explicitly defined in [START_REF] Eliazar | Growth-collapse and decay-surge evolutions, and geometric Langevin equations[END_REF] and z (u) can be obtained from the Lagrange inversion formula as follows: develop u (z) as a power series in z := (z -α) /α, where

u ( z) = z qα + z (qα + pB (α + α z)) maps [0, 1] to [0, 1].
It can in principle be obtained by Faa di Bruno formula for the composition of analytic functions ( [START_REF] Comtet | Analyse combinatoire[END_REF], p. 146), observing z -1 u ( z) = u 1 (u 2 ( z)) where u 1 (z) = 1/ (qα + z) and u 2 ( z) = z (qα + pB (α + α z)) both with well-defined series expansion near 0. Then, by using Lagrange inversion formula ( [START_REF] Comtet | Analyse combinatoire[END_REF], p. 159) we can get z (u) satisfying u ( z (u)) = u and then also the output

E (u τx 0 ,0 ) = 1 α z (u) x0-1 (z (u) -α) = z (u) (α z (u) + α) x0-1 .
The probability generating function z (u) is well-defined in the range u ∈ [0, 1] in the recurrent case and u ∈ [0, 1] → z (u) ∈ [α, 1] is its range.

Harmonic (scale) function and the height of an excursion (recurrent case)

Suppose {X n } , with transition matrix P , is recurrent. Let then P → P with P (0, y) = δ 0,y , forcing state 0 to be absorbing, corresponding to the substitution of processes: X n → Y n := X n∧τx,0 (the process X n started at x 0 and stopped when it first hits 0). The matrix P is stochastic but non irreducible, having {0} as an absorbing class. For the process Y n with transition matrix P , eventually Y ∞ = 0 with probability 1. The harmonic (or scale) sequence φ solves:

P φ = φ, (18)
where φ = (φ (0) , φ (1) , ...) is the right column eigenvector of P associated to the eigenvalue 1. With h 1 and x ∈ {1, ..., h -1} , let

τ x,h = inf (n ≥ 0 : X n ≥ h | X 0 = x) if such an n exists, = +∞ if not.
By induction, with φ (0) := 0 and τ

(h) x := τ x,0 ∧ τ x,h ∀x ∈ {1, ..., h -1} , ∀n ≥ 0: φ (x) = Eφ Y n∧τ x,h = Eφ X n∧τ (h) x .
The harmonic function φ on {1, ..., h}, makes φ X n∧τ (h)

x a martingale. Therefore, as n → ∞, ∀x ∈ {1, ..., h -1}:

φ (x) = Eφ X τ (h) x = φ (0) P (τ x,0 < τ x,h ) + Eφ X τ x,h P (τ x,h < τ x,0 ) . (19)
As h → ∞, both τ x,h and the overshoot X τ x,h → ∞. Assuming there is a solution φ (x) ≥ 0 such that φ (0) = 0 and φ (x) → ∞ as x → ∞, owing to X τ x,h → ∞, yields

P (τ x,h < τ x,0 ) = φ (x) Eφ X τ x,h → h→∞ 0 = P (τ x,∞ < τ x,0 ) , (20)
indeed consistent with the guess: Y n does not escape, it goes extinct with probability 1, whatever the initial condition x is.

Let us now consider the problem of evaluating the law of the height H of an excursion in the recurrent case. Firstly, H = 0 with probability q (corresponding to a trivial excursion). Now, for non-trivial excursions, the event H ≥ h ≥ 1 is realized whenever a first birth event occurs with size β 1 ≥ h or, if 0 < β 1 < h, when the time (starting from β 1 ) needed to first hit {h, h + 1, ...} is less than the time needed to first hit 0, namely the event τ β1,h < τ β1,0 . Hence,

P (H ≥ h) = pP (β 1 ≥ h) + p h-1 x=1 P 0 (β 1 = x) P (τ x,h < τ x,0 ) , where P (β 1 = x) = b x , P (β 1 ≥ h) = y≥h b y and P (τ x,h < τ x,0
) is given in [START_REF] Bingham | Stochastic Processes: Theory and Methods[END_REF]. Note P (H ≥ 1) = p and, as required since P (τ x,∞ < τ x,0 ) = 0, H < ∞ with probability 1.

The transient case

In the transient case (ζ > 1 or ρ > ρ c ), the denominator of Φ x0 (u (z) , z) tends to 0 as u → u (z) and u (z) does not cancel the numerator: u (z) is a true pole of Φ x0 (u, z). When z ∈ [α, 1], u (z) is concave but it has a maximum u (z * ) strictly larger than 1, attained at some z * inside (α, 1), owing to u (α) = 0, u (1) = 1 and u (1) < 0. We anticipate that for some constant C > 0 depending on x 0 ,

P x0 (X n = 0) ∼ C • u (z * ) -n as n → ∞,
stating that {X n } only visits 0 a finite number of times before drifting to ∞.

Large deviations

Consider v (z) := 1/u (z) where u (z), as a pole, cancels the denominator of Φ x0 (u, z) without cancelling its numerator, so

v (z) = qαz + p (z -α) B (z) z -α .
Over the domain 1 ≥ z > α, v (z) is convex with

v (z) = pB (z) -qαα (z -α) -2 and v (z) = pB (z) + 2qαα (z -α) -3 > 0.
Thus v (1) = 1 and v (1) = pαρ -qα α = pρ -qα/α.

We have v (1) > 0 if and only if the chain is transient. In this transient case,

Φ n (z) 1/n → v (z) as n → ∞. Define F (λ) := -log v e -λ = log u e -λ .
The function F (λ) is concave on its definition domain λ ∈ [0, -log α). It is first increasing, attains a maximum and then decreases to -∞ while crossing zero in between. There exists λ * ∈ (0, -log α) such that ω * = F (λ * ) = 0.

With ω ∈ (F (-log α) , F (0)], define

f (ω) = inf 0≤λ<-log α (ωλ -F (λ)) ≤ 0, (21)
the Legendre conjugate of F (λ). The variable ω is Legendre conjugate to λ with ω = F (λ) and λ = f (ω). Note F (-log α) = -∞ and F (0) = v (1) > 0. On its definition domain, f (ω) ≤ 0 is increasing and concave, starting from f (-∞) = -∞ and ending with f (F (0)) = 0 where f (F (0)) = 0. From [START_REF] Varadhan | Large deviations (Special invited paper)[END_REF], we get:

Proposition 7 For those ω in the range [ω * = 0, F (0)] and for any x 0 > 0,

1 n log P x0 1 n X n ≤ ω → n→∞ f (ω) . (22)
In particular, at ω = F (0) = v (1) > 0 where f (F (0)) = -F (0) = 0, we get 1 n X n a.s.

→ v (1) as n → ∞.

To keep ω = F (λ) in the non-negative range [ω * = 0, F (0)], the range of λ should then equivalently be restricted to [0, λ * ]. We clearly have f (ω * ) = -F (λ * ) < 0 and, from

,

- 1 n log P x0 1 n X n ≤ 0 = - 1 n log P x0 (X n ≤ 0) → F (λ * ) as n → ∞. (23)
This shows the rate at which P x0 (X n = 0) decays exponentially with n. Equivalently, with z * = e -λ * , P (X n = 0) ∼ C•u (z * ) -n as guessed and the Green series Φ x0 (1, 0) = G x0,0 (1) = n≥1 P x0 (X n = 0) now is summable for all x 0 ≥ 0 (translating that {X n } visits 0 only a finite number of times). In the transient case,

P (τ x0,0 < ∞) = φ x0,0 (1) = G x0,0 (1)
G 0,0 (1) < 1.

The scale function and the extinction probability

In the transient case, the chain {X n } started at x > 0 can drift to ∞ before it first hits 0. There is thus only a probability smaller than 1 that {X n } gets extinct for the first time.

In the transient case, let then P → P with P (0, y) = δ 0,y forcing state 0 to be absorbing, corresponding to: X n → Y n := X n∧τx,0 . The matrix P is stochastic but non irreducible, having an absorbing class. The harmonic (or scale) sequence φ solves:

P φ = φ, (24)
where φ = (φ (0) , φ (1) , ...) is a column vector. With h 1 and x ∈ {1, ..., h -1} , let

τ x,h = inf (n ≥ 0 : X n ≥ h | X 0 = x) if such an n exists, = +∞ if not.
By induction, with φ (0) = 1 and τ

(h) x := τ x,0 ∧ τ x,h ∀x ∈ {1, ..., h -1} , ∀n ≥ 0: φ (x) = Eφ Y n∧τ x,h = Eφ X n∧τ (h) x .
The harmonic function on {1, ..., h}, makes φ X n∧τ (h) x a martingale, [START_REF] Norris | Markov chains[END_REF]. As n → ∞, ∀x ∈ {1, ..., h -1}:

φ (x) = Eφ (X τx ) = φ (0) P (τ x,0 < τ x,h ) + Eφ X τ x,h P (τ x,h < τ x,0 ) . (25)
As h → ∞, both τ x,h and the overshoot X τ x,h → ∞. Assuming there is a solution φ (x) > 0 such that φ (x) → 0 as x → ∞, yields

P (τ x,0 < τ x,∞ ) = φ (x) , (26)
indeed consistent with the guess: φ > 0 and vanishing at ∞. This expression also shows that φ (x) must be decreasing with x. Note P (τ x,0 < τ x,∞ ) = P (X τx = 0) where τ x := τ x,0 ∧ τ x,∞ . We obtained:

Proposition 8 In the transient case, with φ defined in [START_REF] Harris | The theory of branching processes[END_REF] obeying φ (0) = 1,

φ (x) = P (τ x,0 < τ x,∞ ) ,
is the probability of extinction starting from state x.

Remark: The function φ := 1 -φ clearly also is a (increasing) harmonic function so that so is any convex combination of φ and φ, [START_REF] Neveu | Chaines de Markov et théorie du potentiel[END_REF]. Clearly, φ (x) = P (τ x,∞ < τ x,0 ) = P (X τx = ∞) .

The height of an excursion (transient case)

The law of the height H of an excursion is given by

P (H ≥ h) = P 0 (X 1 ≥ h) + h-1 x=1 P 0 (X 1 = x) P (X τx ≥ h) = P 0 (X 1 ≥ h) + h-1 x=1
P 0 (X 1 = x) P (τ x,h < τ x,0 ) , where P 0 (X 1 = x) = P (0, x) = pb x , P 0 (X 1 ≥ h) = p y≥h b y and, from (25), P (τ x,h < τ x,0 ) = 1 -φ (x) 1 -Eφ X τ x,h .

Note, as required, that H = 0 with probability q and H = ∞ with probability x≥1 P 0 (X 1 = x) φ (x) .

Doob transform: conditioning on non-extinction

In the transient case when the chain can either drift to infinity or go extinct, the harmonic sequence φ plays some additional role in a conditioning. With D φ :=diag(φ (0) , φ (1) , ...), the matrix P φ := D -1 φ P D φ is a stochastic matrix. We have P φ (m, 0) = δ m,0 so that state 0 is inaccessible from any other state than 0 itself: state 0 is isolated and disconnected and so may be removed from the state space. We have ( [START_REF] Rogers | Markov Processes and Martingales[END_REF], page 327, [START_REF] Norris | Markov chains[END_REF]) :

Proposition 9 In the transient case, P φ is the one-step transition matrix corresponding to the process X n∧τx conditioned not to be absorbed at 0. Remark: this is a selection of paths procedure allowing to focus only on those paths of the transient chain that do not go extinct.

Concluding remarks

We revisited the geometric catastrophe model in discrete-time, as a Markovian population dynamics on the non-negative integers. For this process, a collapse move from some state is geometrically distributed so long as it does not exhaust the current value of this state. It is balanced by random growth moves with arbitrary distribution. This process may be viewed as a generalized birth and death chain. Using generating function techniques, conditions under which the birth and death competing events yield a process which is stable (recurrent) have been highlighted. When recurrent, the shape of the invariant probability measure was described. When unstable (transient), the chain either drifts to infinity or go extinct, a feature similar to supercritical branching processes, [START_REF] Harris | The theory of branching processes[END_REF]. The height and length of excursions, extinction probability, time to extinction, have been studied both in the recurrent and transient setups. The harmonic (scale) function was shown to play an important role in the analysis.

Other interesting random walks in the same spirit, but with different collapse rules, were introduced in [START_REF] Neuts | An interesting random walk on the non-negative integers[END_REF]. We plan, as we do here for the truncated geometric case, to lift the veil on some of their intrinsic statistical properties which are expected to be of a completely different nature.

4. 1

 1 No non trivial ( = 0) invariant measure in the transient case Before turning to this question, let us observe the following: suppose ζ > 1 (ρ > ρ c ), so that the super-critical geometric chain is transient. For the same reason as for the null-recurrent case, the chain has no non trivial ( = 0) invariant positive measure either. It is not Harris-transient,[START_REF] Harris | Transient Markov chains with stationary measures[END_REF].

A geometric(α) rv with success probability α takes values in N. A shifted geometric(α) rv with success probability α takes values in N 0 . It is obtained while shifting the former one by one unit.

In the sequel, a boldface variable, say x, will represent a column-vector so that its transpose, say x , will be a row-vector.

Acknowledgments: T.H. acknowledges partial support from the labex MME-DII (Modèles Mathématiques et Économiques de la Dynamique, de l' Incertitude et des Interactions), ANR11-LBX-0023-01. This work also benefited from the support of the Chair "Modélisation Mathématique et Biodiversité" of Veolia-Ecole Polytechnique-MNHN-Fondation X.


Ivan Castillo Camacho

Kai Wang
email: kai.wang@gipsa-lab.grenoble-inp.fr

Data-Dependent Scaling of CNN's First Layer for Improved Image Manipulation Detection

Keywords: Image forensics, convolutional neural network, first-layer convolutional filter, image manipulation detection, stability of variance

Convolutional Neural Networks (CNNs) have become an effective tool to detect image manipulation operations, e.g., noise addition, median filtering and JPEG compression. In this paper, we propose a simple and practical method for adjusting the CNN's first layer, based on a proper scaling of first-layer filters with a data-dependent approach. The key idea is to keep the stability of the variance of data flow in a CNN. We also present studies on the output variance for convolutional filter, which are the basis of our proposed scaling. The proposed method can cope well with different first-layer initialization algorithms and different CNN architectures. The experiments are performed with two challenging forensic problems, i.e., a multi-class classification problem of a group of manipulation operations and a binary detection problem of JPEG compression with high quality factor, both on relatively small image patches. Experimental results show the utility of our method with a noticeable and consistent performance improvement after scaling.

Introduction

Rapid technology development of cameras to capture digital images comes together with a big suite of software to modify an image. Now changes on an image can be so subtle that noticing them for the naked eye is a difficult task. At the same time, the development of techniques that analyze intrinsic fingerprints in image data, i.e., the image forensics research, is one of the most effective ways to solve challenges related to the authenticity of images. The utility of these techniques has taken more importance in recent years because tampered images can now have serious consequences on our modern society, e.g., misleading the public opinion with fake images and presenting falsified images as proof in court.

Basic image manipulation operations such as median filtering, resampling and noise addition are commonly used during the creation of tampered images. Our objective is to detect traces left by these manipulation operations in an image. It is not a surprise that the current trend is to use Convolutional Neural Network (CNN) to detect image manipulations because of its very good forensic performance. In the meanwhile, image forensics researchers in general agree that specific design is required in order to build successful CNNs for forensic tasks. In particular, the CNN's first layer needs to be carefully designed so as to extract useful information relevant to the forensic task at hand [START_REF] Chen | Median filtering forensics based on convolutional neural networks[END_REF][START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF][START_REF] Li | Localization of deep inpainting using high-pass fully convolutional network[END_REF]. Such relevant information is often believed to be in the so-called image residuals, roughly speaking, in the high-frequency components of the image.

Accordingly, for image manipulation detection, a common choice for the CNN's first layer is high-pass filters. Although satisfying results can be achieved, one important aspect, i.e., the stability of the amplitude of the data flow in CNN, has often been ignored or has not been carefully studied. We have the intuition that after the image data passes through a first layer of high-pass filters, the filters' output becomes a weak signal. This would be detrimental to the training of CNN because the data flow shrinks. In this paper, we show that this signal shrinking indeed exists for first-layer filters generated by several popular initialization algorithms that have been used for detecting image manipulations. In addition, with a proper formulation of the first-layer's convolution operation and based on natural image statistics, we provide an intuitive explanation regarding the signal shrinking and subsequently propose a simple scaling method to enhance the output signal. Experimental results, with different first-layer initializations, CNN architectures and classification problems, show that the proposed scaling method leads to consistent improvement of forensic performance.

The remainder of this paper is organized as follows. In Section 2, we briefly review the related work, focusing on the first-layer initialization of CNNs used to detect image manipulations. In Section 3 we present an experimental study, as well as an intuitive theoretical explanation, regarding the variance of the output signal of first-layer filters of different initializations. Our proposed datadependent scaling method is described in Section 4. Experimental results are presented in Section 5. Finally, conclusions are drawn in Section 6.

Related Work

Early methods for detecting basic image manipulation operations were based on feature extraction and classifier training. Different handcrafted features were proposed to detect specific and targeted manipulation operation, e.g., median filtering, resampling and JPEG compression. Afterwards, researchers focus on the more challenging problem of developing a universal method for image manipulation detection. Various methods have been proposed, based on steganalysis features, image statistical model and more recently deep learning.

During the last decade, deep learning methods, including CNNs, have gained outstanding success in a wide range of research problems in the computer vision field. CNNs can learn themselves useful features from given data, effectively replacing the difficult task of handcrafted feature design for human experts. In the recent years, CNNs have also been used to solve image forensic problems. Researchers have noticed a fundamental difference between computer vision tasks and image forensics tasks. The former focuses on the semantic content of images, while the latter often looks for a weak signal representing the difference between authentic and manipulated images. Accordingly, image forensics researchers found that directly applying CNN initialization borrowed from the computer vision field, e.g., the popular Xavier initialization [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF], results in rather limited performance in detecting image manipulation operations [START_REF] Chen | Median filtering forensics based on convolutional neural networks[END_REF][START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF]. Different customized CNN initialization algorithms have been proposed to cope better with forensic tasks. The basic idea of these methods are more or less similar, i.e., generating or using a kind of high-pass filters at the CNN's first layer.

SRM (Spatial Rich Model) filters are one popular and effective way to initialize the first layer of CNNs that are used to solve image forensic problems, e.g., the detection of manipulation operations [START_REF] Chen | Image processing operations identification via convolutional neural network[END_REF], of splicing and copy-move forgeries [START_REF] Liu | Image forgery localization based on multiscale convolutional neural networks[END_REF], and of inpainted images [START_REF] Li | Localization of deep inpainting using high-pass fully convolutional network[END_REF]. SRM filters, a group of handcrafted high-pass filters originally designed for steganalysis [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF], are put at CNN's first layer as initialization and this often leads to very good forensic performances. Indeed, as shown later in this paper, in many cases that we tested, SRM filters outperform other kinds of first-layer filters, especially after the proposed scaling.

Bayar and Stamm [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] proposed a new type of constrained filters for the first layer of a CNN designed to detect image manipulation operations. The idea is to constrain the network's first layer to learn a group of high-pass filters. This is realized by normalizing the filters before each forward pass of the CNN training. The normalization consists of two steps: firstly, the center element of filter is reset as -1; secondly, all the non-center elements are scaled so that they sum up to 1. In this way the sum of all filter elements is 0, and the constrained firstlayer filter behaves like a high-pass one which is effective in suppressing image content. Recently, Castillo Camacho and Wang [START_REF] Castillo Camacho | A simple and effective initialization of CNN for forensics of image processing operations[END_REF] proposed an alternative way of initializing CNN's first layer for image manipulation detection. This is essentially an adaptation of the conventional Xavier initialization [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF] to the situations where it is required to generate high-pass filters after initialization. This method can generate a set of random high-pass filters to be put at CNN's first layer.

In this paper, we consider four different algorithms for first-layer initialization of CNNs with the application to image manipulation detection: the conventional Xavier initialization from the computer vision community [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF], the initialization with SRM filters [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF], Bayar and Stamm's constrained filters [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF], and Castillo Camacho and Wang's high-pass filter initialization [START_REF] Castillo Camacho | A simple and effective initialization of CNN for forensics of image processing operations[END_REF]. We show that all the four methods can produce filters which shrink the input signal at their output, and that our proposed data-dependent scaling can noticeably and consistently improve the forensic performance for all the four initialization algorithms when tested on different CNN architectures and forensic problems.

Variance of Output of Convolutional Filter

It is demonstrated that the stability of the data flow in CNN, as reflected by the variance of the signal in and out a layer, is beneficial for the training of CNN [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF][START_REF] Li | CS231n: Convolutional Neural Networks for Visual Recognition[END_REF]. Ideally, the variance of input and output of a layer should be equal to each other. In this section, we first show that we can predict the variance of the output of a convolutional filter by using statistics of input signal and elements of the filter. Then, we present observations and understandings regarding the output variance for the four initialization algorithms of convolutional filter which are mentioned in the last section. For the sake of brevity, the four algorithms are hereafter called Xavier [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF], SRM [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF], Bayar [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF], and Castillo [START_REF] Castillo Camacho | A simple and effective initialization of CNN for forensics of image processing operations[END_REF].

Formulation

Motivation. We observe potential limitations of the four initialization algorithms and have the intuition that the variance of output of a convolutional filter initialized by them may change substantially compared to the input. SRM [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF] and Bayar [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] do not take into account the output variance during the initialization, because the two algorithms put directly third-party SRM filters or normalized high-pass filters at first layer without modelling the relation between input and output. Xavier [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF] and Castillo [START_REF] Castillo Camacho | A simple and effective initialization of CNN for forensics of image processing operations[END_REF] consider the input-output relation and generate pseudo-random filters. These two initialization algorithms are based on a statistical point of view and realized by drawing pseudo-random samples, so in practice properties of initialized filters may differ for different realizations. Therefore, it is interesting and important to experimentally and theoretically study the actual output variance for each realization of initialized filter.

Formulation for computing output variance. A convolutional layer used in CNN contains a set of learnable filters (also called kernels) [START_REF] Li | CS231n: Convolutional Neural Networks for Visual Recognition[END_REF]. During the forward pass, the kernel moves in a sliding-window manner across the input and computes a weighted sum of the local input data and the kernel. This procedure results in a so-called activation map comprising all the local results computed at every sliding movement of the kernel. Now assume that the kernel contains N scalars denoted by W = (w 1 , w 2 , ..., w N ), then the local input data involved in the computation also contains N scalars, denoted by X = (x 1 , x 2 , ..., x N ). It is easy to see that the local output y is simply the dot product of W and X, as:

y = W, X = N i=1 w i .x i . (1)
In Xavier [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF] and Castillo [START_REF] Castillo Camacho | A simple and effective initialization of CNN for forensics of image processing operations[END_REF], both w i and x i are assumed as independent random variables. In this paper, we take a new and more practical point of view. Since we focus on a proper scaling of a given kernel, we assume that the kernel elements w i are known constants, which can be generated by any initialization algorithm. In addition, we do not consider x i as independent; instead, we consider them as mutually correlated random variables reflecting the natural image statistics [START_REF] Simoncelli | Natural image statistics and neural representation[END_REF]. With these assumptions and based on the property of variance of weighted sum of variables, we can compute the variance of the output y as

Var(y) = Var N i=1 w i .x i = N i=1 N j=1 w i w j Cov(x i , x j ) = N i=1 w 2 i Var(x i ) + 2 1≤i <j≤N w i w j Cov(x i , x j ).
(

) 2
Table 1. Considered manipulation operations and their parameters.

Median filtering

F ilterSize = 3 Gaussian blurring StandardDeviation = 0.5, F ilterSize = 3 Additive Gaussian noise StandardDeviation = 1.1 Resampling
ScalingF actor ∈ {0.9, 1.1}

JPEG compression

QualityF actor ∈ {90, 91, ..., 100}

The last expression just divides all the relevant terms into two groups: variance terms and covariance terms of the input signal components (x 1 , x 2 , ..., x N ). Furthermore, it is well-known that natural images have approximate translation invariance [START_REF] Simoncelli | Natural image statistics and neural representation[END_REF], implying that Var(x i ), i = 1, 2, ..., N are almost identical. In addition, the neighboring pixels are usually highly-correlated [START_REF] Simoncelli | Natural image statistics and neural representation[END_REF], which means that Cov(x i , x j ) is close to Var(x i ). We approximate Var(x i ) by Var(x), the overall variance of input. Then we have the following approximation of Eq. ( 2):

Var(y) ≈ Var(x)   N i=1 w 2 i + 2 1≤i <j≤N w i w j C ij   , (3)
with C ij = Cov(x i , x j )/Var(x) which are in practice smaller than but very close to 1 for small natural image patches due to high correlation of neighboring pixels.

Experimentally the above equation approximates very well the output variance.

It also helps us to understand the output variance of popular initialized filters used for manipulation detection, as presented in the remaining of this section.

Convolutional filter initialized with high-pass filter

We first consider convolutional filter initialized as each of the 30 SRM filters1 of shape 5×5 (so here N = 25). In order to test on real data for convolutional filter, we take as input 64 × 64 grayscale image patches generated from the Dresden database [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF]. The image manipulation operations that we want to detect are listed in Table 1. We then compute the variance of output of each SRM filter by two different methods: the first one with Eq. ( 3) and the second one with actual convolution between the input and the filter. Hereafter, we call the first as covariance-based method because Eq. ( 3) is mainly based on the covariance terms Cov(x i , x j ) of the input signal components (x 1 , x 2 , ..., x N ), and we call the second one as convolution-based method. For the first method, the covariance terms are estimated from 5 × 5 small patches (same size as SRM filters) which are randomly extracted from the aforementioned 64×64 Dresden image patches. The results of Var(y)/Var(x), i.e., the ratio of output and input variance, are shown in Fig. 1. We can see that the amplitude of Var(y)/Var(x) is very small for all 30 SRM filters, which lies basically in a range from 0 to 0.016 with a mean of about 0.005. The output of majority of SRM filters has a variance smaller than 1% of input variance, reflecting the signal shrinkage. It can also be observed that the two methods to obtain output variance give very close results of Var(y)/Var(x), implying the coherence of the prediction by Eq. ( 3) with the practical convolution results.

In order to understand the small output variance for SRM filters, we start from one important property of these high-pass filters. Like many high-pass filters, e.g., Laplacian filter, the sum of all filter elements is equal to 0 for all 30 SRM filters of shape 5×5 (i.e., N = 25), which means that we have N i=1 w i = 0 (cf., link in footnote 1). It is then easy to deduce that N j=1 w j . N i=1 w i = N i=1 N j=1 w i .w j = 0. By dividing the w i .w j terms into two groups, we obtain

N i=1 w 2 i + 2 1≤i <j≤N w i w j = 0. (4)
The left-hand side of the above equation is almost same as the term in the bracket of Eq. ( 3), except that in the above Eq. ( 4) we replace C ij by 1. As mentioned earlier, for small natural image patches, we have the property that C ij are usually smaller than but very close to 1. This is verified by experiments on Dresden database where the minimum C ij value is 0.9573 for 5 × 5 small patches. Not surprisingly, this minimum C ij value is attained between two pixel positions which are the farthest from each other within the 5 × 5 small patch.

From the above analysis we can see that the term in the bracket of Eq. ( 3) is close to 0, which results in a small value of output variance Var(y) for 30 SRM filters. This intuitively explains the small Var(y)/Var(x) values shown in Fig. 1.

Regarding the high-pass filter initialized by Bayar [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] and Castillo [START_REF] Castillo Camacho | A simple and effective initialization of CNN for forensics of image processing operations[END_REF], we simulated 10, 000 filters with both algorithms and calculated the variance of output of simulated filters. We also observe very small values of the ratio of output-input variance, with 0.005 and 0.006 as mean value of Var(y)/Var(x), respectively for Bayar and Castillo. Due to space limit, we do not show detailed results, but these small mean values further confirm the behaviour that high-pass filters result in small output variance with significant signal shrinkage.

Convolutional filter with Xavier initialization

With curiosity, we also carry out studies for Xavier [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF] which generates 5×5 filters filled with pseudo-random samples drawn from a zero-mean uniform distribution. We created 10, 000 Xavier filters using PyTorch and Fig. 2 shows two histograms of occurrences of output-input variance ratio, i.e., Var(y)/Var(x): the left one is for the range of 0 to 3 with a bin width of 0.02, while the right one shows detailed occurrences for the first bin of the left histogram for the range of 0 to 0.02 with a bin width of 0.001. We computed the mean value of Var(y)/Var(x) for the 10, 000 simulations of Xavier and found that the mean is close to 1 (desired value of Xavier); this is because of a long tail of big values that we do not completely show in Fig. 2. The left histogram of Fig. 2 does not have a peak around 1; instead, the highest occurrences occur near 0. This is a little surprising yet understandable according to Eq. (3). In fact, the elements of Xavier filter are drawn from a zero-mean distribution, so the bracket term in Eq. ( 3) tends to have a small value. However, due to numerical sampling and in particular considering the relatively low number of 25 pseudo-random samples (for 5 × 5 filter), it is possible for the bracket term to take big values in certain simulations. Experimentally this bracket term of Eq. ( 3) can be as big as 13 for some Xavier filters. In addition, from the right histogram of Fig. 2, for Xavier the occurrences of Var(y)/Var(x) being very small values, i.e., less than 0.01, is very low: 109 out of 10, 000 simulations (i.e., around 1% probability). In contrast, the majority of this variance ratio is less than 0.01 for high-pass filters as presented in last subsection. We guess that it is still related to the numerical sampling of Xavier because it can be rare to have 25 pseudo-random samples which sum up to a value extremely close to 0.

Scaling of Convolutional Filter

From results and analysis in Section 3, we can see that the variance of output of convolutional filter initialized by popular algorithms can be significantly smaller than the variance of input. This is particularly true for high-pass filter: the ratio of output-input variance Var(y)/Var(x) is usually smaller than 0.01. The output signal after convolution operation substantially shrinks. This can be detrimental to the training of CNN, and as shown later in Section 5 the CNN training sometimes fails in such situations. Using a data-dependent approach (i.e., dependent on input data), we propose a simple yet effective scaling of the first-layer convolutional filter. The idea is to keep the variance stable after scaling for the input and output of any given filter generated by popular initialization algorithms. Corresponding to the two methods to compute output variance in Section 3, we propose two different ways to calculate the scaling factor s, as presented below. After obtaining the scaling factor, the elements of the given filter W = (w 1 , w 2 , ..., w N ) are properly scaled as W = s.W. We then initialize the first-layer filter with the scaled version W. Covariance-based method. From Eq. ( 3), it can be seen that in order to make Var(y) and Var(x) approximately equal to each other, we need to compensate for the effect of the term in the bracket. So the scaling factor is computed as:

s = 1   N i=1 w 2 i + 2 1≤i <j≤N w i w j C ij   . (5)
In practice, we take random small patches of the same shape of the convolutional filter to be scaled (e.g., 5 × 5) from a small portion of the training data. We then estimate the variance and covariance terms on these small patches to obtain the values of C ij = Cov(x i , x j )/Var(x). Afterwards the scaling factor s is computed by using Eq. ( 5), and at last we obtain the scaled version W of any given filter W from the considered four initialization algorithms.

Convolution-based method. This is a straightforward approach. The output y is computed, for a small portion of the training data x as input, by carrying out the convolution operation. The scaling factor is simply calculated as

s = Var( x)/Var( y). (6)
From a practical point of view, the covariance-based method might be a slightly better option than the convolution-based method mainly because of its higher flexibility. In fact, for the covariance-based method, the computation of the variance and covariance terms of the input can be performed only once for any number of filters for which we want to scale. By contrast, the convolutionbased method has to be rerun every time we have a new filter to analyze. Nevertheless, it is worth mentioning that both methods are experimentally quick enough to be used in CNN initialization. The running time is about several seconds, as presented in the next paragraph.

According to our experiments, for a training set of about 100, 000 images of 64×64 pixels from all classes, taking 10% of the training data for the convolutionbase method and 10 small patches (e.g., of 5 × 5 pixels) per image of the 10% training data for the covariance-based method, we achieve a good trade-off between computation time and stability of the result. Using more training data has very small impact on the obtained scaling factor. Even using 100% of the training set results in a change smaller than 0.1%. The amount of time to calculate the scaling factor is less than 3 seconds per first-layer filter for both methods, on a desktop with Intel R Xeon E5-2640 CPU and Nvidia R 1080 Ti GPU (covariancebased method on CPU and convolution-based method on GPU). This is run for one time before the CNN training. The computation time increases very slowly when having more filters for the covariance-based method, because as mentioned above the variance and covariance terms can be reused. We believe that the computation time of scaling factor is negligible when compared to the typical time required to train a CNN model.

Experimental Results

Several experiments are performed in order to test and show the efficiency of our proposed scaling. These experiments consider the four filter initialization algorithms mentioned earlier, two CNN architectures (CNN of Bayar and Stamm [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] and a smaller CNN without fully-connected layer designed by ourselves), and two forensic problems (a multi-class problem of detecting a group of manipulation operations and a binary problem of detecting JPEG compression of high quality factor). For the multi-class problem, we also consider a different number of filters used in the first layer of the CNN of Bayar and Stamm [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF]. The implementation and experiments were conducted using PyTorch v1.4.0 with Nvidia R 1080 Ti GPU. The experimental data was created from the Dresden database [START_REF] Gloe | The Dresden image database for benchmarking digital image forensics[END_REF]. Full-resolution Dresden images are split for training, validation and testing with ratio of 3:1:1 and converted to grayscale. Patches of 64 × 64 pixels were randomly extracted from full-resolution grayscale Dresden images. This relatively small size of image patches makes the forensic problems more challenging.

Multi-class problem with CNN of Bayar and Stamm [1]

We first consider the multi-class problem of classifying six different kinds of image patches: the original patches and the five classes of manipulated patches as explained in Table 1. The parameters for the resampling and JPEG compression manipulations are taken randomly from the specified sets in Table 1. The total number of patches in training set is 100, 000 (≈ 16, 667 patches per class), while the number of patches in testing set is 32, 000 (≈ 5, 333 patches per class). The number of training and testing samples is same as in [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF]. It is worth mentioning that the manipulations and their parameters in this paper are borrowed from [START_REF] Castillo Camacho | A simple and effective initialization of CNN for forensics of image processing operations[END_REF] and more challenging than those in [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF]. The patch size is also smaller than [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF]: our patches are of 64 × 64 pixels, while [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] mainly considers 256 × 256 patches. We use the successful CNN architecture of Bayar and Stamm [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] in this set of experiments and initialize the three filters in the CNN's first layer with four different algorithms: Bayar [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF], SRM [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF], Castillo [START_REF] Castillo Camacho | A simple and effective initialization of CNN for forensics of image processing operations[END_REF], and Xavier [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. We carry out 5 runs for each algorithm and the corresponding two scaled versions. For SRM, for each run we randomly select 3 filters from the pool of 30 SRM filters. We compare each original initialization algorithm with their scaled versions obtained with the covariance-based method and the convolution-based method presented in Section 4. For fair comparisons, we make sure that for each run the scaled versions share the same "base filters" of the original version before performing scaling. We follow exactly the same training procedure described in [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF], including number of epochs, optimization algorithm, learning rate schedule, etc.

For Bayar algorithm [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF], we have tested two variants of the scaling of the firstlayer filters. The first one ("Bayar A.") follows closely the idea of Bayar's original constrained training strategy: we carry out scaling of the normalized high-pass filter at the beginning of each forward pass (please refer to the second last paragraph of Section 2 for detail of the normalization procedure proposed in [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF]). The second variant ("Bayar B.") is computationally cheaper and less complex: the scaling of normalized high-pass filter is only performed in the initialization, and we no longer impose normalization constraint during training. Our intuition behind the second variant is that with a proper scaling of initialized filters even a free training without the constraint of [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] may provide satisfying performance.

The detection performances in terms of test accuracy (i.e., classification accuracy on testing set) for this multi-class problem are presented in Table 2. The reported results are average of 5 runs with randomness, e.g., different first-layer "base filters". However, for each run, the "base filters" are the same for the original and scaled versions: original version direct uses these filters, while scaled versions apply proper scaling on the "base filters" and then use the scaled ones.

We can see from Table 2 that the test accuracy of all the initialization algorithms is consistently and noticeably improved after scaling. Improvement of at least 1.13% and as high as 2.74% is obtained. We also observe that the results of the two scaling methods are very close to each other. We checked the computed scaling factors and found that they are indeed almost identical for the two methods. Furthermore, for the two scaling variants of Bayar [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF], vari- ant B gives slightly better results, which is also computationally cheaper as it only performs scaling in initialization without enforcing any constraint during CNN training. This implies that with a good initialization after proper scaling, it might not be necessary to impose training constraint. It is worth mentioning that the results of Bayar in Table 2 are in general lower than those reported in [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] because we now consider a more challenging forensic problem with more difficult manipulations and on smaller patches. The results of Castillo in Table 2 are better than those presented in [START_REF] Castillo Camacho | A simple and effective initialization of CNN for forensics of image processing operations[END_REF]. This may be due to the differences in the number of training epochs (we follow [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] and train with more epochs) and in the adopted optimization algorithm and learning schedule. In addition, the three kinds of high-pass filters (especially SRM) indeed outperform Xavier, before and after scaling. This demonstrates the difference between forensics and computer vision tasks. Nevertheless, the performance of Xavier is also improved after scaling because as analyzed in Section 3.3 Xavier can also result in small variance of output. In fact, for Xavier the probability to have Var(y) smaller than half of Var(x) is about 52.20% in our 10, 000 simulations.

We also observe that the proposed scaling helps to have quicker increase of forensic performance during CNN training. We show in Fig. 3 curves of test accuracy of SRM and Xavier (average of 5 runs), before and after the covariancebased scaling. It can be observed that the convergence speed is considerably improved for SRM. For both algorithms, the curve of scaled version is always above that of original version during the whole 60 epochs. It is also interesting to notice that the scaled Xavier performs slightly better than the original SRM.

With 30 filters at first layer. Next, we present results for the same multi-class problem while changing the number of filters in the first layer of the CNN of [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] to 30. We make this change for two reasons: first, to test our approach with a different number of filters in the first layer; and second, to use all the 30 SRM filters which is a common practice in image forensics, e.g., for detecting splicing Table 3. Test accuracy for the multi-class forensic problem (in %, average of 5 runs). We still use the CNN of [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] but change the number of first-layer filters from 3 to 30. and copy-move forgeries [START_REF] Liu | Image forgery localization based on multiscale convolutional neural networks[END_REF]. For this scenario we still test the four initialization algorithms but only use variant B for scaled Bayar as it proved to obtain slightly better results while being computationally cheaper. The results are presented in Table 3. Again, we observe that scaling the filters with any of the two methods leads to consistently better test accuracy, with an improvement ranging from 1.13% to 4.23%. We notice from Tables 2 and3 that after increasing the number of first-layer filters, 1) the original version of high-pass initialization (Bayar, Castillo and SRM) has slightly improved or comparable performance while the accuracy of Xavier decreases; and 2) the scaled version of Bayar, Castillo and SRM has comparable performance with the case of 3 filters while Xavier has noticeable improvement. We guess the reason for the good performance of scaled Xavier may be that with 30 filters there is more chance to have a very good filter which after scaling can improve the result. Understanding these observations is not the focus of our paper, and we plan to conduct further analysis in the future.

Initialization

JPEG binary problem with CNN of Bayar and Stamm [1]

We notice in the multi-class problem that JPEG compression is the most difficult manipulation to detect. In this section we consider the binary classification between original patches and JPEG compressed patches with parameters in Table 1 (i.e., very high quality factor between 90 and 100). This allows us to test the proposed scaling on a different challenging forensic problem. We use the CNN of Bayar and Stamm [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF]. The number of training and testing patches per class is the same as in last subsection. All CNN training settings are kept unchanged.

For this binary problem, we consider two initialization algorithms of Bayar [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] and SRM [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF], original and scaled versions (variant B for scaled Bayar). Table 4 presents the obtained results (average of 5 runs). This challenging problem makes the original version of both Bayar and SRM struggle to achieve a good performance. Especially, training of SRM can occasionally fail, leading to accuracy close to random guess. Much better average test accuracy is achieved by scaled versions. For SRM [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF], a boost of more than 14% is obtained with scaling.

We show in Fig. 4 some curves of test accuracy for scaled (covariance-based) and original SRM [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF]. The curve of scaled version shows the maximum and minimum test accuracy together with the average at each epoch among the 5 runs. For the original version of SRM we can have very different results. Therefore, we show the best and the worst curves of test accuracy among all the 5 runs. As we can see the worst case does not improve during the whole procedure and the test accuracy remains close to 50%. The difference may come from the randomly selected three first-layer SRM filters in each run (certain SRM filters perform worse than others according to our observation). We would like to mention that for each run, although we select randomly different SRM filters, the same selected filters are used to carry out comparisons between the original and scaled versions. Therefore, even for filters that result in bad performance for the original version, we can obtain a much better performance after scaling them.

Multi-class and binary problems on a different smaller CNN

We then test both the multi-class and JPEG binary problems on a different CNN designed by ourselves. We first describe the architecture of this smaller network. Let Ck(M or A) denote a Convolutional-BatchNorm-Tanh(-MaxPool or -AveragePool) layer with k filters. For the first layer we use Hk which denotes a Convolutional layer with k filters. The architecture of our smaller CNN is H3-C40M-C25M-C20M-C15M-C6A. The first four layers have a kernel size of 5 × 5 while for the last two layers the kernel size is 1×1. All convolutional stride size is 1. The first layer and the last two layers do not have zero-padding, for the other layers the padding size is 2. This is a network without fully-connected layer. To compare with, the architecture proposed by Bayar and Stamm [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] is H3-C96M-C64M-C64M-C128A-F200-F200-F6, where Fk denotes a fully-connected layer with k neurons and Tanh. The number of learnable parameters of the CNN of [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] is about 337K, while our smaller CNN has about 41K parameters. Using our smaller CNN, we test both the multi-class and the binary problems on a different CNN architecture. All the data preparation and experimental setting are the same as those described in Sections 5.1 and 5.2. For this set of experiments, we compare the original and scaled versions of Bayar [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF] and SRM [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF] (variant B for scaled Bayar). Table 5 presents the obtained results. We can see that in all cases the scaled version leads to improved performance compared to the original version. The improvement of test accuracy goes from 0.93% for multi-class problem with Bayar, to 5.27% for binary problem with SRM.

Our objective in this subsection is to show that with a different CNN, our proposed scaling can still reliably improve the performance for different initialization algorithms and forensic problems. Meanwhile, it can be noticed that performance is better with our smaller CNN when compared to the network of [START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF]. The understanding of this point is beyond the scope of this paper. Our guess is that the forensic problems and/or the amount of data cope better with the smaller CNN's size (less parameters) and architecture (only comprising convolutional layers without fully-connected layer). The thorough analysis and understanding of the relationships between these factors is one part of our future work.

Conclusion

We propose a new and effective scaling approach for adjusting first-layer filters of CNNs used for image manipulation detection. The proposed scaling is computationally efficient and data-dependent (i.e., scaling factor dependent on the input). We also present theoretical and experimental studies which help to understand why the ratio of output-input variance for first-layer convolutional filter can be a (very) small value. Experimental results, with different CNNs, filter initialization algorithms and forensic problems, show that our proposed scaling can consistently improve performance of CNN-based image manipulation detection. Although practically and intuitively we can understand the effectiveness of the proposed filter scaling operation, a rigorous theoretical analysis would be necessary to explain the observed performance improvement. We would like to conduct research on a possible scaling of deeper layers and to extend experiments to more CNNs (e.g., XceptionNet and ResNet), training settings and multimedia security problems. It is also interesting to study the impact of amount of training data and CNN architecture on the forensic performance. One limitation of our work is that we only consider image manipulation detection under ideal laboratory conditions. We plan to carry out relevant studies for more challenging and practical forensic problems, e.g., the detection of malicious and realistic image forgeries and image forensics with mismatch between training and testing data.

Fig. 1 .

 1 Fig. 1. The value of Var(y)/Var(x) for each of the 30 SRM filters obtained by using the covariance-based method of Eq. (3) and the convolution-based method.

Fig. 3 .

 3 Fig. 3. Curves of test accuracy (average of 5 runs) for the multi-class forensic problem, during the whole 60 training epochs of the CNN of [1]. The curves are for SRM and Xavier, original version and scaled version by the covariance-based method.

Fig. 4 .

 4 Fig. 4. Curves of test accuracy for the JPEG binary forensic problem: scaled version of SRM (average of 5 runs) with bars of maximum and minimum accuracy at each epoch among 5 runs; and the best run and the worst run of original version of SRM.

  Two histograms of occurrences of output-input variance ratio Var(y)/Var(x) for 10, 000 Xavier filters. Please refer to main text for detailed explanation.

	600
										40
	Count 400									20 30
	200
										10
	0									0
	0.0	0.3	0.6	0.9	1.2	1.5 Var(y)/Var(x) 1.8	2.1	2.4	2.7	3.0	0.000	0.002	0.004	0.006	0.008 Var(y)/Var(x) 0.010 0.012	0.014	0.016	0.018	0.020
	Fig. 2.

Table 2 .

 2 Test accuracy for the multi-class forensic problem (in %, average of 5 runs). The experiments were performed with four initialization algorithms and their scaled versions for first-layer filters of the CNN of Bayar and Stamm[START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF]. In parentheses is the improvement of scaled version compared to the corresponding original version.

	Initialization	Original version	Convolution-based scaling	Covariance-based scaling
	Bayar [1] A. Bayar [1] B.	94.19	96.04 (+1.85) 96.15 (+1.96)	96.02 (+1.83) 96.22 (+2.03)
	Castillo [2]	93.71	96.45 (+2.74)	96.42 (+2.71)
	SRM [5]	94.39	96.54 (+2.15)	96.55 (+2.16)
	Xavier [7]	93.48	94.61 (+1.13)	94.71 (+1.23)

Table 4 .

 4 Test accuracy for the binary JPEG forensic problem (in %, average of 5 runs). The experiments were performed with Bayar and SRM, original and scaled versions (variant B for scaled Bayar), on the CNN of[START_REF] Bayar | Constrained convolutional neural networks: A new approach towards general purpose image manipulation detection[END_REF].

	Initialization	Original version	Convolution-based scaling	Covariance-based scaling
	Bayar [1] B.	88.27	90.80 (+2.53)	90.80 (+2.53)
	SRM [5]	78.24	92.33 (+14.09)	92.44 (+14.20)

Table 5 .

 5 Test accuracy for the multi-class and binary problems with our proposed smaller CNN without fully-connected layer (in %, average of 5 runs). The columns of "Bayar" and "SRM" present results of original version. "Scaling-conv" and "Scalingcov" represent respectively convolution-based and covariance-based scaling method.

	Problem Bayar Scaling-conv Scaling-cov SRM Scaling-conv Scaling-cov
	Multi-class 95.24 96.17(+0.93) 96.18(+0.94) 95.72 97.06(+1.34) 97.09(+1.37)
	Binary	90.56 93.72(+3.16) 93.74(+3.18) 89.85 95.11(+5.26) 95.12(+5.27)

The 30 SRM filters can be found in the class of SrmFiller, starting from line 347 of this webpage https://github.com/tansq/WISERNet/blob/master/filler.hpp

Acknowledgments

This work is partially funded by the French National Research Agency (DEFALS ANR-16-DEFA-0003, ANR-15-IDEX-02) and the Mexican National Council of Science and Technology (CONACYT).