Lionel Ronsoux
email: lionel.ronsoux@insa-lyon.fr

Mariette Moevus
email: mariette.moevus@gmail.com

Yves Jorand
email: yves.jorand@insa-lyon.fr

Sandrine Maximilien
email: sandrine.maximilien@insa-lyon.fr

Christian Olagnon
email: christian.olagnon@insa-lyon.fr

Romain Anger
email: romain.anger@grenoble.archi.fr

Laetitia Fontaine
email: laetitia.fontaine@grenoble.archi.fr

Poured Earth as concrete

Keywords: poured earth, clay concrete, grain packing, clay dispersion

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Earthen construction provides a number of advantages: local raw materials, low energy manufacturing, interesting hygrothermal behaviour. However, for a more extensive use, it can sometimes be necessary to optimize the manufacturing process or to improve the mechanical behaviour of this material as well as its durability.

For 2000 years, earthen construction techniques have hardly changed. However, at the scientific level, the various constituents of earth (aggregates, clay and water) have provided some of the most innovative subjects of scientific investigation of the last twenty years. Compared with dry and wet granular media [START_REF] Bocquet | Moisture induced ageing in granular media and the kinetics of capillary condensation[END_REF], earth brought an additional level of physicochemical complexity at the finest scale of clays [START_REF] Van Damme | L'eau et sa représentation[END_REF]. The emergence of nanoscience and the physical chemistry of colloids and interfaces will undoubtedly contribute to the understanding of phenomena occurring in these materials: what is the nature of the forces responsible for the cohesion of particles? What is the role of surface interactions between the clay particles in the viscoelastic and rheological behaviour of the saturated material? What is the role of capillary water in the cohesion of earth?

These scientific questions arise for other materials such as, for example, cementitious materials [START_REF] Pellenq | Why does concrete set? The nature of cohesion Forces in hardened cement-based materials[END_REF]. Indeed clays present many similarities with the calcium silicate hydrate, the main cement hydrate [START_REF] Viallis-Terrisse | Interaction des Silicates de Calcium Hydratés, principaux constituants du ciment, avec les chlorures d'alcalins[END_REF]. From the technological and industrial points of view, outstanding advances also prepare a favorable ground for a more innovative use of earth materials. For instance, the enhanced properties of self-compacting and self-leveling concretes are partly due to the use of cement dispersing agents and a better distribution of the granular packing [START_REF] Vernet | Ultra-durable concretes: structure at the micro-and nanoscale[END_REF][START_REF] Flatt | The Rheology of Cementitious Materials[END_REF]. Similarly, in the field of industrial ceramics, effective dispersants and binders are used for shaping before firing.

The improvement of the properties of earthen materials can benefit from expanded research and findings applicable to industrialised materials. Following this idea, we applied the two main advances cited before (optimisation of granular packing and use of dispersing agents) to earthen materials.

The purpose of this paper is to present the first experimental results obtained through this research program. The effect of dispersing agents on the rheological behaviour of soil suspensions is presented. Then mortars are prepared and characterised by uniaxial compression tests. The results show that the use of dispersing agents leads to higher fluidity and compressive strength.

MATERIALS AND EXPERIMENTS

Soil

Natural soil from a quarry in Dauphiné (France) has been used in this study. This soil is suitable for rammed earth construction. It was sieved at 100µm. The particle size distribution is presented in Figure 1. A mineral analysis by X-ray diffraction mentions the presence of the following elements: quartz, clays (illite and muscovite) and alkali compounds (albite and microcline). The specific density has been measured with a helium pycnometer: 2.61 g/cm 3 . The soil is stocked at ambient temperature and ambient hygrometry in the laboratory.

Standardized siliceous sand, with a specific density of 2.65 g/cm 3 , is used to prepare mortars and compressive samples. Its particle size distribution is presented in Figure 1.

Additives

Two polymer additives have been selected in this paper for their dispersing effect: (1) Na-HMP: sodium hexametaphosphate Na 6 O 18 P 6 , and (2) PAA: sodium salt of polyacrylic acid (C3H4O2) n with a molecular weight of 2100 g/mol. Both of them are widely used in the ceramic industry for their dispersing and deflocculating properties.

The dispersing mechanisms that can be expected by using such additives are the following: adsorption of the molecules at the surface of clay particles, electrostatic repulsion between the particles with adsorbed molecules, steric repulsion between the chains of macromolecules adsorbed on clays (Lewis, 2000, Fontaine andAnger, 2009). These electrostatic and steric forces keep the particles well dispersed.

Dispersant adsorption is highly dependent on the electrostatic interactions between the dispersant and the clay particle surface. Thus the surface chemistry of the solid phase and the solution properties of the polyelectrolyte are important parameters, and these parameters depend on the pH and the ionic strength of the solution. These considerations oblige to work properly in well-defined laboratory conditions, using pure water with controlled quality, in order to observe reproducible results.

Rheological characterization of dilute suspensions

Some of the authors of this paper have shown interest in the rheological approach to understand the physico-chemical mechanisms involved in the implementation of earth in a plastic state [START_REF] Anger | Influence de la teneur en sel et du pH sur la plasticité du matériau terre[END_REF]. The same approach is used here to characterize the effect of dispersants on the rheological behaviour of earth.

As it is expected that these polymers will interact only with the finest particles, we chose to work on suspensions containing only the soil fraction smaller than 100µm. It is also expected that the action of the polymers will be exacerbated in a dilute suspension rather than in a concentrated paste (mud). Thus, we decided to study the relationship between the rheological behaviour and the polymer content on dilute suspensions containing a 36% vol. of solid particles in water.

We start by preparing a large batch of dilute suspension of earth. Some reverse osmosis water is stirred in a beaker with a magnetic shaker. The earth is sifted and slowly poured into the vortex of water. When all the earth has been added, the beaker is covered with a plastic film to limit evaporation. The suspension is stirred for 24 hours.

After 24 hours, several 80 mL samples are prepared with different concentrations of dispersant: 0%, 0.05%, 0.2%, 0.5%, 0.8%. These mixtures are stirred again for 5 minutes for homogenization, and then characterized by rheology.

Rheological characterization is carried out with a co-axial viscosimeter Haake VT500. The liquid suspension is introduced between the two coaxial cylinders of the viscosimeter. The external cylinder is fixed and the internal cylinder can rotate either at a constant rotation speed (Ω) or at constant shear stress (τ). The shear rate γ' (or gradient of speed inside the fluid) is directly proportional to Ω.

The shear rate is increased from 0 to 1000 s -1 within 2.5 min, maintained constant during 30 s, then decreased to 0 s -1 in 2.5 min. We obtain a rheogram, in other words a diagram showing the shear stress versus the shear rate of the suspension analyzed. The studied supensions behave like Bingham plastic liquids, and their rheogram can be described by the Bingham relation which is a two-parameter model :

τ = τ s + η plast .γ'
where τ s is the yield stress below which the fluid cannot be sheared, and η plast is the plastic viscosity : the lower it is, the easier is to shear the liquid, once the threshold is exceeded. The weaker these parameters are, the more fluid the suspension is.

The rheological behaviour is determined on the basis of the descending rheogram. Such a procedure allows the samples to be characterized in the same agitated state. Two measures are performed for each sample to verify that the measure is reproducible and that the suspension is not altered during the experiment.

Optimisation of particle size distribution

Earthen materials used for construction are commonly described according to their granulometry: an earth may be qualified as sandy, fine, rocky, silty, clayey, etc. Each earth is a mix of grains with different sizes and the proportions change from one earth to another.

The packing of solid particles is of essential importance to the understanding of granular materials used in different size classes and related problems which appear in many fields of science and industrial processes, such as in civil engineering or ceramics.

A proper packing of grains is really important, first of all to improve the placement of the material, reduce shrinkage, facilitate the drying process and secondly to improve its final properties like its compressive strength or abrasion resistance. The optimization of particle size distribution (PSD) is the way to achieve a good packing.

In cement concrete, ingredients with a wide range of PSDs are combined. Therefore, different approaches can be found for the composition of concrete mixtures. Many authors have worked on this topic since the 1900s and have offered various solutions. The mix design concept discussed here and applied to earth consists in the proportioning of different grain sizes to formulate a continuously graded material. This idea is realized by the formulation of a mix using the equation of cumulative volume fraction from Andreasen: CPFT = (D/D max ) q

(1)

with CPFT being the "Cumulative Percent Finer Than", D being the particle size, D max being the maximum particle size and q being the distribution moduli.

Particle size distributions for the fine fraction of Dauphiné soil, the siliceous sand and the mortar mix for mechanical characterization are summarized in the Figure 1. The mortar mix elaborated to produce samples for compressive tests is based on soil fines and standardized siliceous sand. In order to obtain a very good particle packing, the proportions of each constituent are defined in accordance with Andreasen's model. This model, one of the most commonly used in the field of ceramics, indicates a distribution moduli between 1/3 and 1/2 to reach the best packing [START_REF] Dinger | One dimensional packing of spheres, Part II[END_REF]. Our mix shows a good fitting with the distribution corresponding to q=0.37, noted "Model" in Figure 1. This mix corresponds to 70% in weight siliceous sand plus 30% soil fines. In this particular case, proportions in weight are equivalent to proportions in volume due to the close specific densities of raw materials.

Uniaxial compressive tests

Conventional mechanical tests concerning building materials remain simple and usually measure only one characteristic. The uniaxial compressive test is very common, but the majority of compressive strength values given in the scentific literature are inaccurate because the sample is not loaded in purely uniaxial compression [START_REF] Fontaine | Cohésion et comportement mécanique de la terre comme matériau de construction[END_REF][START_REF] Morel | Compressive strength testing of compressed earth blocks[END_REF][START_REF] Mollion | Etude du comportement mécanique du pisé[END_REF]. If the compressive test is perfectly carried out, with no restraint from the machine plates and with direct strain measurement during the test, the intrinsic behaviour of the material can be assessed accurately.

To achieve this, we must overcome the edge effects induced by the friction between the sample and the compression plates. It is also necessary that the surfaces of the sample be parallel and smooth. Therefore, in this study, great care has been taken in the preparation of samples and the execution of compression tests.

Three different mortars are prepared with the soil/sand mix, with the optimised particle size distribution presented before: (1) reference mortar without dispersing agent, (2) with PAA, (3) with Na-HMP. The polymer content introduced in mortars ( 2) and ( 3) is equal to 0.5% in weight compared to the mass of fine soil. This quantity has been chosen on the basis of the rheological results presented below.

The fine particles are first mixed with freshly filtered reverse osmosis water and the dispersant. The preparation is mixed for 1 hour in a mortar mixer, then packed in a sealed bag. It is then kept for 24 h at 20°C. Then this mixture is again introduced into the mixer and the standard sand is incorporated. The mortar is mixed for 5 min at low speed and 5 min at high speed. The amount of water is adjusted to obtain a similar consistency for the three mortars. This consistency corresponds to a spread of 15.5 cm after 15 shocks on a standard shock table.

Cylindrical samples (diameter 30mm, height 55mm) are prepared by casting the preparation with shocks to help the material to spread. The samples are demoulded after 24 hours and left to dry at room temperature in the laboratory.

A fine earthen mortar is applied after one week on the surfaces of the cylinders to obtain very smooth and parallel surfaces. The samples are placed in a climatic chamber (20°C and 50% relative humidity) until cons tant mass is reached, typically within a few days.

The uniaxial compression tests are performed on an INSTRON 8562 electromechanical machine. A roller is placed above the specimen to make up for parallelism defects. Oiled plastic wrap is placed between the specimen and each compression steel plate, in order to reduce friction. A LVDT sensor +/-10mm measures the displacement between the compression plates. An extensometer (gauge length 25mm) is placed directly in contact with the specimen and measures the deformation in the central portion of the sample, where the stress field is homogeneous.

RESULTS

Rheology

Viscosity measurement results for Dauphiné soil suspensions are presented below on Fig. 2 and Fig. 3. The effects of Na-HMP and PAA on the viscosity of the suspensions are very similar. The optimum polymer content is 0.5% for both additives. The plastic visosity is divided by 3 with both dispersing agents, which represents an important reduction of viscosity. The yield stress is drastically reduced with small amounts of additives, about 0.5% in weight. This drastic reduction of yield stress indicates that the suspensions which exhibit a shear thinning behaviour without a dispersing agent, operate a huge change and exhibit, at the optimum of the dispersing agent, a behaviour close to a newtonian fluid. These important reductions of both rheological parameters at the same water content result in a real difference in the fluidity of the suspensions. This remarkable similitude between the effects of both additives on the earthen suspension studied here suggests that the mechanisms of dispersion involved are very similar. This was not expected at such a high level due to the complex nature and numerous different constituants in the fine fraction of Dauphiné soil.

These rheological results are of great interest for the future. They show the potential of possibilities to transfer the physico-chemical knowledge from other fields involved with the very fine particles.

Mechanics

The three mortars presented above were characterized in uniaxial compression. shows that the curves obtained for the four samples from the same series are very reproducible. This helps to validate the quality of the tests done.

The results from the compression tests performed on the three materials are given in Tab.1, and the Fig. 5 shows the effect of the two dispersants used on the water content w init , the dry mass ρ dry , and the compressive strength f c . In this figure, the values are normalized to facilitate comparisons. We see that the addition of PAA or of Na-HMP allows reducing the amount of water necessary for the implementation of earth mortars, at a defined consistency. It also shows that the dry density is slightly increased by the addition of the dispersant. The most striking result is the strong increase in terms of compressive strength: it is multiplied by 1.8 for the mortar containing PAA, and by 1.6 for the mortar containing Na-HMP, as compared to the reference mortar.

Such an increase in resistance cannot be due only to the very small increase in density. It is very likely that the use of dispersants leads to a different microstructure, with a more homogeneous distribution of pores and a better organization of clay platelets, for example. Microstructural characterizations are being made to better understand the origin of this phenomenon.

These first results seem very promising. They open an interesting path for the implementation of unstabilized earth concretes, without the need for a hydraulic binder.

CONCLUSION

The goal of the study presented herein is to transfer some of innovations from the fields related to cementitious materials and industrial ceramics and apply them to earthen construction. First, the authors optimized the granular packing of their soil so as to obtain a solid material with low porosity. Then, they characterized the influence of dispersants on the soil's rheological and mechanical properties.

Experimental results show on one hand that the use of very small amounts of polymer dispersants can greatly reduce the viscosity and yield stress of dilute suspensions of soil. In the case of an earth mortar with a granular optimization, this results in a lower amount of water needed for its implementation, to obtain the desired consistency. Thus, the addition of a dispersant could help reduce the drying time for an earth concrete. On the other hand, the mechanical characterization of earth mortars shows that the addition of a dispersant can greatly increase the compressive strength of earth. This result is particularly interesting because it opens new ways to explore and achieve unstabilized earth concretes, with a resistance that could match the mechanical properties of earth concrete stabilized with a hydraulic binder, without affecting the hydro-regulating properties of natural earth.

  Figure 1 Particle size distributions

Figure 2

 2 Figure 2 Plastic viscosity of suspensions with 36% vol. fraction of soil.

Figure 3

 3 Figure 3 Threshold shear stress value of suspensions with 36% vol. fraction of solid particles.

  Fig 4.


Hanany Ould-Baba
email: hanany.ould-baba@utc.fr

Vincent Robin
email: vincent.robin@utc.fr

Jérôme Antoni
email: jerome.antoni@insa-lyon.fr

Concise formulae for the cumulant matrices of a random vector

Keywords: Moments, cumulants, cumulant matrix, commutation matrix, independent component analysis, rank of cumulant matrices, moment matrix of a Gaussian vector

Concise formulae are given for the cumulant matrices of a random vector up to order 6. In addition to usual matrix operations, they involve only the Kronecker product, the vec operator, and the commutation matrix. Orders 5 and 6 are provided here for the first time; the same method as provided in the paper can be applied to compute higher orders. An immediate consequence of these formulae is to return 1) upper bounds on the rank of the cumulant matrices and 2) the expression of the sixth-order moment matrix of a Gaussian vector. Due to their conciseness, the proposed formulae also have a computational advantage as compared to the repeated use of Leonov and Shiryaev formula.

Introduction

Cumulants are useful quantities for characterizing some important statistical properties, such as independence, gaussianity, or possible (a)symmetry of probability distributions. They were first introduced in 1899 by T.N. Thiele under the name halfinvariants (see Hald [START_REF] Hald | The early history of the cumulants and the gram-charlier series[END_REF] where original Thiele's paper is translated into English). In particular, higher order cumulants make it possible to solve problems that are otherwise intractable on the second order. This is typically the case in Independent Component Analysis (ICA, as introduced in Comon [START_REF] Comon | Independent component analysis, a new concept?[END_REF]), where fourth-order cumulants are used to blindly separate (over-)determined (as many as or less sources than sensors) mixtures of independent sources (e.g. Cardoso and Soulomiac [START_REF] Cardoso | Blind beamforming for non gaussian signals[END_REF], Albera et al. [START_REF] Albera | Icar: a tool for blind source separation using fourth-order statistics only[END_REF]. Recently, several researchers have investigated the use of sixth-order cumulants (Albera et al. [START_REF] Albera | Sixth order blind identification of underdetermined mixtures (birth) of sources[END_REF]) and, more generally, cumulants of order 2k, k > 2, (Albera et al. [START_REF] Albera | Blind identification of over-complete mixtures of sources (biome)[END_REF]) for the separation of under-determined mixtures of sources (more sources than sensors). This paper is concerned with the matrix expression of the cumulants of a random vector. Loosely speaking, the cumulant matrix of a random vector is a collection of the various crosscumulants between its elements. It is often convenient to arrange such a collection into a matrix (a square symmetric matrix in the case of even order cumulants) for both notational and practical reasons. This makes it easier to to handle cumulants for algorithmic perspectives and to exploit some of their statistical properties, such as "matrix redundancy" (see e.g. Albera et al. [START_REF] Albera | Sixth order blind identification of underdetermined mixtures (birth) of sources[END_REF][START_REF] Albera | Icar: a tool for blind source separation using fourth-order statistics only[END_REF][START_REF] Albera | Blind identification of over-complete mixtures of sources (biome)[END_REF]). An important example is ICA where the blind separation of statistically independent sources is solved by zeroing all cross-cumulants between independent sources: if properly set up, the issue then boils down to the diagonalization of a cumulant matrix (see Cardoso and Soulomiac [START_REF] Cardoso | Blind beamforming for non gaussian signals[END_REF], Kollo [START_REF] Kollo | Multivariate skewness and kurtosis measures with an application in ica[END_REF]).

It is a well-known fact that scalar cumulants can be evaluated from lower-order moments by use of the classical Leonov and Shiryaev formula (see Leonov and Shiryaev [8] or Speed [START_REF] Speed | Cumulants and partition lattices1[END_REF], Mac Cullagh [START_REF] Cullagh | Tensor Methods in Statistics[END_REF]). This trivially applies to all (cross-)cumulants of a vector which can eventually be rearranged into a matrix. The main drawback of that approach is to involve lengthy formulae and, when calculation is of concern, a computational burden, which both become rapidly prohibitive as the dimension of the random vector increases. It is therefore advisable to search for direct and more concise expressions of cumulant matrices. One solution is actually returned by Neudecker matrix derivative (Neudecker [11]) (see e.g. Kollo and von Rosen [12, pages 187-188]). Unfortunately, it faces two difficulties which limit its practical use:

(i) the resulting cumulant matrix is not square except on the second order;

(ii) matrix derivatives of the logarithm of the characteristic function are hardly tractable for orders higher than four and, according to the authors knowledge, no formulae have ever been derived for those cases.

Limitation (i) can be overcome on the fourth order by using an algebraic technique such as proposed by Jinadasa and Tracy [START_REF] Jinadasa | Higher order moments of random vectors using matrix derivatives[END_REF] and Tracy and Sultan [START_REF] Tracy | Higher order moments of multivariate normal distribution using matrix derivatives[END_REF] to obtain the matrix expression of the fourth-order moment of a Gaussian vector; the same result is found in Kollo [START_REF] Kollo | Multivariate skewness and kurtosis measures with an application in ica[END_REF], Loperfido [START_REF] Loperfido | A note on the fourth cumulant of a finite mixture distribution[END_REF]. However, point (ii) has not been solved till now. In this paper we give the expressions of fifth and sixth-order cumulant matrices of a real-valued random vector. The expression of the sixthorder cumulant matrix is symmetric.

The paper is organized as follows. Section 2 resumes the basic definitions of moments and cumulants of a random vector and deduces Leonov and Shiryaev formula as a consequence of de Faà di Bruno formula (following the lines in Lukacs [START_REF] Lukacs | Applications of faa di bruno's formula in mathematical statistics[END_REF] and Hardy [START_REF] Hardy | Combinatorics of partial derivatives[END_REF]). The main result of the paper -a theorem providing the formulae of cumulant matrices up to order six -is given in section 3. At the same time a general strategy is described as how to obtain similar results to any order. In section 4, the formulae are used to get, without pain, i) upper bounds on the rank of cumulant matrices and ii) the sixth-order moment matrix of a Gaussian vector. Since the latter is a known formula (see for instance Magnus and Neudecker [START_REF] Magnus | The commutation matrix: Some properties and applications[END_REF], Tracy and Sultan [START_REF] Tracy | Higher order moments of multivariate normal distribution using matrix derivatives[END_REF], and more recently Schott [START_REF] Schott | Kronecker product permutation matrices and their application to moment matrices of the normal distribution[END_REF]), it somehow validates the expression of the cumulant matrix from which it is derived.

Finally, section 5 concludes on the advantages implied by the proposed matrix expressions as compared to the classical use of Leonov and Shiryaev formula.

Generalities on moment and cumulant matrices

This section first recalls the link between the scalar moments and cumulants of a random vector, in order to pave the way for our main theorem given in section 3.

Scalar moments and cumulants

Let x = (x 1 , . . . , x p ) ′ be a real-valued random vector of dimension p × 1. Its characteristic function φ is defined on R p as

φ(u) = E exp(i(u, x)) = E         p j=1 exp(iu j x j )         . (1)
Let n ∈ N; vector x is said to admit moments of order n if

E [|x| n ] < ∞, (2)
where |x| = p j=1 x 2 j 1/2 stands for the euclidian norm of R p . In that case, (i) φ is n times differentiable and (ii) all the scalar (cross-)moments of order k ≤ n of x are well defined as

(∀( j 1 , . . . , j k ) ∈ [[1, p]] k ) mom (k) (x j 1 , . . . , x j k ) def = E x j 1 • • • x j k = (-i) k ∂ k φ ∂u j 1 • • • ∂u j k (0). (3)
Moreover, since φ is a continuous function on R p and φ(0) = 1, it is necessarily non-zero in a neighborhood of 0 ∈ R p . The second characteristic function (also referred to as the "cumulant function" in the literature),

ψ(u) = log (φ(u)) , (4)
is thus considered. By composition under condition (2), it is n times differentiable in a neighborhood of 0; this makes possible the definition of scalar (cross-)cumulants of order k ≤ n of x as

(∀( j 1 , . . . , j k ) ∈ [[1, p]] k ) cum (k) (x j 1 , . . . , x j k ) déf = (-i) k ∂ k ψ ∂u j 1 • • • ∂u j k (0). (5)
Since Eq. ( 4) implicitly connects moments to cumulants, it suggests two methods to compute the later in terms of the former:

1. expand up to order n the following function

log          1 + n k=1 (i) k k! ( j 1 ,..., j k )∈[[1,p]] k mom (k) (x j 1 , . . . , x j k )u j 1 • • • u j k          (6)
where the sum in the brackets is the Taylor series to order n of φ(u);

2. use de Faà di Bruno formula which returns the partial derivatives of a composition of applications.

The first method is carried out for instance in Mac Cullagh [START_REF] Cullagh | Tensor Methods in Statistics[END_REF]. However, calculation becomes rapidly intractable as the order increases: Mac Cullagh provides formulae only up to order six; a few authors, like Albera and Comon [START_REF] Albera | Asymptotic performance of contrast-based blind source separation algorithms[END_REF], have attempted to pursue the calculations up to order eight. The second method has been used by Lukacs [START_REF] Lukacs | Applications of faa di bruno's formula in mathematical statistics[END_REF] in the case of random variables (p = 1); using a generalization of de Faà di Bruno formula to the multivariable case, several authors have proposed general formulae to calculate cumulants of a random vector in terms of its moments, which only differ in the way terms are collected. One of the simplest is provided by Hardy [START_REF] Hardy | Combinatorics of partial derivatives[END_REF]; with the notations of the present paper, it reads

∂ k ∂u j 1 • • • ∂u j k f • φ(u) = π∈P([[1,k]]) f (|π|) • φ(u) B∈π ∂ |B| l∈B ∂u j l φ(u) (7)
where f is a function of one variable, P([ [1, k]]) denotes the set of partitions of [ [1, k]], and |E| the cardinality of set E such that, for instance, with B = {2; 3; 6},

∂ |B| l∈B ∂u j l = ∂ 3 ∂u j 2 ∂u j 3 ∂u j 6 . (8)
The main advantage of Eq. ( 7) is to be fully developed, that is it no longer contains factorials assigned to the collection of equal partial derivatives when some indices j l happen to coincide.

In order to compute cumulants in terms of moments, one sets f = log, u = 0 and φ(0) = 1 in Eq. [START_REF] Kollo | Multivariate skewness and kurtosis measures with an application in ica[END_REF]. In particular,

f (n) (1) = (-1) n-1 (n -1)! (9)
and, in accordance with (3),

∂ |B| l∈B ∂u j l φ(0) = i |B| mom (|B|) (x j l : l ∈ B). (10)
By noting that for any π ∈ P([ [1, k]]) it comes B∈π |B| = k, one eventually arrives at the Leonov and Shiryaev formula:

cum (k) (x j 1 , . . . , x j k ) = π∈P([[1,k]]) (-1) |π|-1 (|π| -1)! B∈π mom (|B|) (x j l : l ∈ B). (11)
Alternatively, in order to compute moments in terms of cumulants, one sets f = exp, ψ in place of Φ, u = 0 and ψ(0) = 0 in Eq. [START_REF] Kollo | Multivariate skewness and kurtosis measures with an application in ica[END_REF], so as to obtain

mom (k) (x j 1 , . . . , x j k ) = π∈P([[1,k]]) B∈π cum (|B|) (x j l : l ∈ B). (12)
This formula will be used in section 4.2 to find the expression of the sixth-order moment, m 6 (x), of a centered Gaussian random vector x.

Moment and cumulant matrices

There are many different ways to arrange the p k scalar moments of Eq.( 3) in a matrix. It is natural to select those arrangements which are maximally symmetric, as is often practiced in the literature (see e.g. Schott [START_REF] Schott | Kronecker product permutation matrices and their application to moment matrices of the normal distribution[END_REF], Jinadasa and Tracy [START_REF] Jinadasa | Higher order moments of random vectors using matrix derivatives[END_REF], Tracy and Sultan [START_REF] Tracy | Higher order moments of multivariate normal distribution using matrix derivatives[END_REF]).

Definition 1 (The moment matrix). The moment matrix of a real-valued random vector x = (x 1 , . . . ,

x p ) ′ , E [|x| n ] < ∞, n ∈ N, of even order 2k ≤ n is defined as m 2k (x) = E ⊗ k xx ′ (dimension p k × p k ) ( 13
)
and that of odd order, 2k + 1 ≤ n, as

m 2k+1 (x) = E ⊗ k xx ′ ⊗ x (dimension p k+1 × p k ), ( 14
)
where ⊗ is the Kronecker product and

⊗ k xx ′ = xx ′ ⊗ • • • ⊗ xx ′ k matrices p×p .
Note there are several ways to write down the moment matrix since, for any pair of vectors a and b of dimension p × 1,

ab ′ = a ⊗ b ′ = b ′ ⊗ a. (15)
It is easy to verify that all scalar (cross-)moments are included in definition 1. To do so, let decompose x = j∈[ [1,p]] x j e j onto the canonical basis of R p ; thus, from the multilinearity property of the Kronecker product and the linearity of the expected value,

m 2k (x) = E ⊗ k xx ′ = E          ⊗ k ( j 1 , j 2 )∈[[1,p]] 2 x j 1 x j 2 e j 1 e ′ j 2          = E          ( j 1 , j 2 )∈[[1,p]] 2 x j 1 x j 2 e j 1 e ′ j 2 ⊗ • • • ⊗ ( j 2k-1 , j 2k )∈[[1,p]] 2 x j 2k-1 x j 2k e j 2k-1 e ′ j 2k          = ( j 1 , j 2 ,..., j 2k-1 , j 2k )∈[[1,p]] 2k E x j 1 x j 2 • • • x j 2k-1 x j 2k e j 1 e ′ j 2 ⊗ • • • ⊗ e j 2k-1 e ′ j 2k = ( j 1 , j 2 ,..., j 2k-1 , j 2k )∈[[1,p]] 2k mom (2k) (x j 1 , . . . , x j 2k )e j 1 e ′ j 2 ⊗ • • • ⊗ e j 2k-1 e ′ j 2k . (16)
Similarly, for odd orders,

m 2k+1 (x) = ( j 1 , j 2 ,..., j 2k-1 , j 2k , j 2k+1 )∈[[1,p]] 2k+1 mom (2k+1) (x j 1 , . . . , x j 2k+1 )e j 1 e ′ j 2 ⊗ • • • ⊗ e j 2k-1 e ′ j 2k ⊗ e j 2k+1 . (17)
It is then possible to define the cumulant matrix in a similar way as the moment matrix.

Definition 2 (The cumulant matrix). Let x = (x 1 , . . . , x p ) ′ be a real-valued random vector of dimension p × 1, such that E [|x| n ] < ∞ for a given n ∈ N. The cumulant matrix of even order, 2k ≤ n, of x is a p k × p k matrix defined as

c 2k (x) = ( j 1 , j 2 ,..., j 2k-1 , j 2k )∈[[1,p]] 2k cum (2k) (x j 1 x j 2 • • • x j 2k-1 x j 2k )e j 1 e ′ j 2 ⊗ • • • ⊗ e j 2k-1 e ′ j 2k (18)
whereas the cumulant matrix of odd order, 2k + 1 ≤ n, is a p k+1 × p k matrix defined as

c 2k+1 (x) = ( j 1 , j 2 ,..., j 2k-1 , j 2k , j 2k+1 )∈[[1,p]] 2k+1 cum (2k+1) (x j 1 x j 2 • • • x j 2k-1 x j 2k x j 2k+1 )e j 1 e ′ j 2 ⊗ • • • ⊗ e j 2k-1 e ′ j 2k ⊗ e j 2k+1 . (19)
Except for the special case of a Gaussian random vector, cumulants are much more difficult to compute than moments. Hitherto, the only way to computing the cumulant matrix c n (x) has consisted in 1. deriving the p n scalar cumulants from the scalar moments of orders lower than n via Eq. ( 11), 2. arranging the scalar cumulants in a matrix according definition 2. However, Eq. ( 11) tends to return lengthy expressions that are hardly tractable beyond order four or so; for instance, the sixth-order cumulant requires the consideration (and computation) of 41 moments. It is the object of this paper to compute the cumulant matrix c n (x) directly from the moment matrices of orders lower than n, thus providing formulae that are much more concise and convenient to handle.

Cumulant matrices from lower-order moment matrices

Statement of results for orders 2 to 6

In order to state the theorem of the paper, it is first necessary to introduce the commutation matrix and the vec operator. The following definitions and properties are borrowed from Magnus and Neudecker [START_REF] Magnus | The commutation matrix: Some properties and applications[END_REF].

Lemma 1 (The commutation matrix). Let p ≥ 1 be an integer. The commutation matrix K pp of dimension p 2 × p 2 is defined as

K pp = (i, j)∈[[1,p]] 2 e i e ′ j ⊗ e j e ′ i , (20)
where (e i ) i∈[ [1,p]] stands for the canonical basis of R p .

Let us define two other related commutation matrices, of dimension p 3 × p 3 :

P 3,2 = I p ⊗ K pp (21) P 2,1 = K pp ⊗ I p . (22)
Therefore, if a, b, and c are three vectors of dimension p × 1, the following equalities hold:

K p (a ⊗ b) = b ⊗ a (23) P 3,2 (a ⊗ b ⊗ c) = a ⊗ c ⊗ b (24) P 2,1 (a ⊗ b ⊗ c) = b ⊗ a ⊗ c. ( 25
)
Since K ′ pp = K pp , it comes P ′ 3,2 = P 3,2 , P ′ 2,1 = P 2,1 , and

(a ′ ⊗ b ′ )K p = b ′ ⊗ a ′ (26) a ′ ⊗ b ′ ⊗ c ′ P 3,2 = a ′ ⊗ c ′ ⊗ b ′ (27) a ′ ⊗ b ′ ⊗ c ′ P 2,1 = b ′ ⊗ a ′ ⊗ c ′ . ( 28
) Proof. K p (a ⊗ b) = (i, j)∈[[1,p]] 2 e i e ′ j ⊗ e j e ′ i (a ⊗ b) = (i, j)∈[[1,p]] 2
(e i e ′ j a)

a j e i ⊗ (e j e ′ i b) b i e j = (i, j)∈[[1,p]] 2 b i a j e i ⊗ e j = b ⊗ a (29) P 3,2 (a ⊗ b ⊗ c) = I p ⊗ K pp (a ⊗ b ⊗ c) = a ⊗ K pp (b ⊗ c) = a ⊗ c ⊗ b (30) P 2,1 (a ⊗ b ⊗ c) = K pp ⊗ I p (a ⊗ b ⊗ c) = K pp (a ⊗ b) ⊗ c = b ⊗ a ⊗ c (31)
We are now in a position to state the main result of the paper.

Theorem 2. Let x = (x 1 , . . . , x p ) ′ be a real-valued random vector of dimension p×1, with assumed zero-mean for simplicity (E [x] = 0). Let us further assume that E [|x| n ] < ∞ for a given integer n ≥ 6 so as to guarantee the existence of moments and cumulants up to order six. Next let us define the four matrices:

K = I p 2 + K pp (32) P = I p 3 + P 2,1 + P 3,2 P 2,1 (33) R = I p 3 + P 3,2 + P 2,1 P 3,2 (34) Q = I p 3 + P 3,2 + P 2,1 + P 2,1 P 3,2 + P 3,2 P 2,1 + P 2,1 P 3,2 P 2,1 (35)
which depend only on dimension p. Thus, the cumulant matrices of orders 2 to 6 read

c 2 (x) = m 2 (x) = E xx ′ (36) c 3 (x) = m 3 (x) = E xx ′ ⊗ x (37) c 4 (x) = m 4 (x) -K m 2 (x) ⊗ m 2 (x) -vecm 2 (x)vec ′ m 2 (x) (38) c 5 (x) = m 5 (x) -R m 3 (x) ⊗ m 2 (x) K -P m 3 (x) ′ ⊗ vecm 2 (x) -vecm 3 (x)vec ′ m 2 (x) (39) c 6 (x) = m 6 (x) -Γ (4,2) (x) -Γ (3,3) (x) + 2Γ (2,2,2) (x), (40) 7
respectively, where

Γ (4,2) (x) = R m 4 (x) ⊗ m 2 (x) R ′ + M 4 (x) ⊗ vec ′ m 2 (x) P ′ + P M 4 (x) ′ ⊗ vecm 2 (x) Γ (3,3) (x) = R m 3 (x) ⊗ m 3 (x) ′ P ′ + vecm 3 (x)vec ′ m 3 (x) Γ (2,2,2) (x) = Q m 2 (x) ⊗ m 2 (x) ⊗ m 2 (x) + P m 2 (x) ⊗ vecm 2 (x)vec ′ m 2 (x) P ′ M 4 (x) = E xx ′ ⊗ x ⊗ x .

Outline of the proof

The sketch of the proof consists of (a) injecting Eq. ( 11) that returns cumulants in terms of moments into Eqs. ( 18) and ( 19) which define cumulant matrices; (b) expressing the terms of the resulting sum with moment matrices.

Let us detail the first step in the case of even orders (the principle is identical for odd orders):

c 2k (x) = ( j 1 ,..., j 2k )∈[[1,p]] 2k         π∈P([[1,2k]]) (-1) |π|-1 (|π| -1)! B∈π mom (|B|) (x j l : l ∈ B)         e j 1 e ′ j 2 ⊗ • • • ⊗ e j 2k-1 e ′ j 2k = s∈P(2k) (-1) |s|-1 (|s| -1)!Γ s (x) (41)
where P(2k) denotes the set of all s = (s 1 , . . . , s l ), where s 1 , . . . , s l are strictly positive integers such that s 1 + . . . + s l = 2k partitions, |s| denotes the length of s -that is the number l of integers in s -and

Γ s (x) = π∈P s ([[1,2k]])          ( j 1 ,..., j 2k )∈[[1,p]] 2k B∈π mom (|B|) (x j l : l ∈ B)e j 1 e ′ j 2 ⊗ • • • ⊗ e j 2k-1 e ′ j 2k          T π , (42)
For the second step, it is convenient to browse set P(2k) in the lexicographic order, as suggested in Eq. ( 43), so that the first term in Eq. (41) will correspond to s = (2k), |s| = 1, P s ([ [1, 2k]

]) = {{{1; • • • ; 2k}}}; that is partition π = {{1; • • • ; 2k}} with one block B = {1; • • • ; 2k}, such that Γ (2k) (x) = T {{1;•••;2k}} = ( j 1 ,..., j 2k )∈[[1,p]] 2k mom (2k) (x j l : l ∈ B)e j 1 e ′ j 2 ⊗ • • • ⊗ e j 2k-1 e ′ j 2k = m 2k (x). ( 45
)
In addition, in each product B∈π of (42), blocks B will be browsed from the smallest to the largest, as suggested by Eq. (43). Therefore, for each partition s of the integer 2k, there exists a unique type-s partition π of set [ [1, 2k]] that globally preserves the natural order on [ [1, 2k]]. This will be denoted as the canonical partition associated with s. For instance, the canonical partition of type s = (2, 2) is π s = {{1; 2}; {3; 4}}. Since the blocks of a canonical partition are consecutive, the corresponding terms are easily found. In general, if the type-s partition is such that the number of odd indices is not greater than the number of even indices plus one, then the corresponding term in ( 42) is simply

T π s = ( j 1 ,..., j 2k )∈[[1,p]] 2k B∈π s mom (|B|) (x j l : l ∈ B)e j 1 ⊗ e ′ j 2 ⊗ • • • ⊗ e j 2k-1 ⊗ e ′ j 2k = B∈π s m |B| (x). ( 46
)
For non-canonical partitions, e j 1 ⊗ e ′ j 2 ⊗ • • • ⊗ e j 2k-1 ⊗ e ′ j 2k will be rearranged according to the blocks in partitions π s . Unfortunately, there is no general matrix formulation in this case. However, since e j ⊗ e ′ k = e ′ k ⊗ e j , we will make use of the identity

e j 1 ⊗e ′ j 2 ⊗• • •⊗e j 2k-1 ⊗e ′ j 2k = e j 1 ⊗• • •⊗e j 2k-1 ⊗e ′ j 2 ⊗• • •⊗e ′ j 2k = (e j 1 ⊗• • •⊗e j 2k-1 )(e j 2 ⊗• • •⊗e j 2k ) ′ . ( 47
)
For orders greater than 4, lemma 1 will be used to rearrange Kronecker products of length 3. For blocks where the number of odd indices is neither equal to nor greater than the number of even indices, the corresponding factor will not be of type m |B| (x). Since orders up to 6 are considered only in this paper, this will be fixed by using the vec operator and the moment matrix M 4 (x) = E [xx ′ ⊗ x ⊗ x], and its transpose. Finally, let us note that for a zero-mean random vector, all terms involving a first-order moment -corresponding to partitions with blocks of unit lengthwill be nil.

Proof for order 3

The only partition of 3 which does not involve 1 is (3):

P (3) ([[1, 3]]) = {{{1; 2; 3}}}. Thus c 3 (x) = Γ (3) (x) = ( j 1 , j 2 , j 3 )∈[[1,p]] 3 mom (3) (x j 1 , x j 2 , x j 3 )e j 1 e ′ j 2 ⊗ e j 3 = m 3 (x). (48)

Proof for order 4

The only partitions of 4 (see Eq. ( 43)) which do not involve 1 are (4) and (2, 2):

P (4) ([[1, 4]]) = {{{1; 2; 3; 4}}} P (2,2) ([[1, 4]]) = {{{1; 2}; {3; 4}}; {{1; 3}; {2; 4}}; {{1; 4}; {2; 3}}} .
In order to avoid any ambiguity, each term of Γ (2,2) (x) will be dentoted by T π where π is the corresponding partition in (42). Therefore,

c 4 (x) = m 4 (x) -Γ (2,2) (x) (49) = m 4 (x) -m 2 (x) ⊗ m 2 (x) T {{1;2};{3;4}} - ( j 1 ,..., j 4 )∈[[1,p]] 4
mom (2) (x j 1 , x j 3 )mom (2) (x j 2 , x j 4 )e j 1 e ′ j 2 ⊗ e j 3 e ′ j 4

T {{1;3};{2;4}} - ( j 1 ,..., j 4 )∈[[1,p]] 4
mom (2) (x j 1 , x j 4 )mom (2) (x j 2 , x j 3 )e j 1 e ′ j 2 ⊗ e j 3 e ′

For e j 1 e ′ j 2 ⊗ e j 3 e ′ j 4 = e j 1 ⊗ e ′ j 2 ⊗ e j 3 ⊗ e ′ j 4 = e j 1 ⊗ e j 3 ⊗ e ′ j 2 ⊗ e ′ j 4 = (e j 1 ⊗ e j 3 ) ⊗ (e j 2 ⊗ e j 4 ) ′ (51) = K pp (e j 3 ⊗ e j 1 )(e j 2 ⊗ e j 4 ) ′ = K pp (e j 3 e ′ j 2 ⊗ e j 1 e ′ j 4 ), (52) it comes

T {{1;3};{2;4}} =          ( j 1 , j 3 )∈[[1,p]] 2 mom (2) (x j 1 , x j 3 )e j 1 ⊗ e j 3                   ( j 2 , j 4 )∈[[1,p]] 2 mom (2) (x j 2 , x j 4 )e j 2 ⊗ e j 4          ′ = vecm 2 (x)vec ′ m 2 (x), (53)
T {{1;4};{2;3}} = K pp          ( j 2 , j 3 )∈[[1,p]] 2 mom (2) (x j 2 , x j 3 )e j 3 e ′ j 2          ⊗          ( j 1 , j 4 )∈[[1,p]] 2 mom (2) (x j 1 , x j 4 )e j 1 e ′ j 4          = K pp m 2 (x) ⊗ m 2 (x). ( 54
)
Upon collecting m 2 (x) ⊗ m 2 (x), one finally arrives at

c 4 (x) = m 4 (x) -Km 2 (x) ⊗ m 2 (x) -vecm 2 (x)vec ′ m 2 (x). ( 55
)
3.5. Proof for order 5 The only partitions of 5 which do not involve 1 are ( 5) et [START_REF] Cardoso | Blind beamforming for non gaussian signals[END_REF][START_REF] Comon | Independent component analysis, a new concept?[END_REF]: One may check that P (3,2) ([ [START_REF] Hald | The early history of the cumulants and the gram-charlier series[END_REF][START_REF] Albera | Sixth order blind identification of underdetermined mixtures (birth) of sources[END_REF]]) contains 5! 3!2! = 10 elements. As before, each term of Γ (3,2) (x) will be denoted as T π where π is the corresponding partition in Eq. (42). Thus

P (5) ([[1, 5]]) = {{{1
c 5 (x) = m 5 (x) -Γ (3,2) (x), ( 56
)
where

Γ (3,2) (x) = π∈P (3,2) ([[1,5]]) T π .
Next, the T π 's are computed by factoring terms of order 3 in first positions. The principle is illustrated for instance on T {{2;3;5};{1;4}} :

T {{2;3;5};{1;4}} = ( j 1 ,..., j 5 )∈[[1,p]] 5
mom (3) (x j 2 , x j 3 , x j 5 )mom (2) (x j 1 , x j 4 )e j 1 ⊗ e ′ j 2 ⊗ e j 3 ⊗ e ′ j 4 ⊗ e j 5 = ( j 1 ,..., j 5 )∈[[1,p]] 5 mom (3) (x j 2 , x j 3 , x j 5 )mom (2) (x j 1 , x j 4 )(e j 1 ⊗ e j 3 ⊗ e j 5 ) ⊗ (e ′ j 2 ⊗ e ′ j 4 )

= ( j 1 ,..., j 5 )∈[[1,p]] 5
mom (3) (x j 2 , x j 3 , x j 5 )mom (2) (x j 1 , x j 4 )P 2,1 P 3,2 (e j 3 ⊗ e j 5 ⊗ e j 1 ) ⊗ (e ′ j 2 ⊗ e ′ j 4 )

= P 2,1 P 3,2 ( j 1 ,..., j 5 )∈[[1,p]] 5 mom (3) (x j 2 , x j 3 , x j 5 ) =mom (3) (x j 3 ,x j 2 ,x j 5 )
mom (2) (x j 1 , x j 4 )(e j 3 ⊗ e ′ j 2 ⊗ e j 5 ) ⊗ (e j 1 ⊗ e ′ j 4 )

= P 2,1 P 3,2 m 3 (x) ⊗ m 2 (x). ( 57
)
Following similar lines, (66)

T {{1;2;3};{4;5}} = m 3 (x) ⊗ m 2 (x) (58) T {{1;2;4};{3;5}} = m 3 (x) ′ ⊗ vecm 2 (x) (59) T {{1;3;4};{2;5}} = m 3 (x) ⊗ m 2 (x)K pp (60) T {{2;3;4};{1;5}} = P 2,1 m 3 (x) ′ ⊗ vecm 2 (x) (61) T {{1;2;5};{3;4}} = P 3,2 m 3 (x) ⊗ m 2 (x) (62
The last step is to collect terms in m 3 (x) ⊗ m 2 (x), m 3 (x) ′ ⊗ vecm 2 (x), and vecm 3 (x)vec ′ m 2 (x), respectively:

Γ (3,2) (x) = I p 3 + P 3,2 + P 2,1 P 3,2 R m 3 (x) ⊗ m 2 (x) + I p 3 + P 3,2 + P 2,1 P 3,2 R m 3 (x) ⊗ m 2 (x)K pp + I p 3 + P 2,1 + P 3,2 P 2,1 P m 3 (x) ′ ⊗ vecm 2 (x) + vecm 3 (x)vec ′ m 2 (x) = Rm 3 (x) ⊗ m 2 (x) I p 2 + K pp K +Pm 3 (x) ′ ⊗ vecm 2 (x) + vecm 3 (x)vec ′ m 2 (x) = Rm 3 (x) ⊗ m 2 (x)K + Pm 3 (x) ′ ⊗ vecm 2 (x) + vecm 3 (x)vec ′ m 2 (x). (67)
3.6. Proof for order 6 The 4 partitions of 6 which do not involve 1 are: (6), (4, 2), (3, 3) et (2, 2, 2). Therefore c 6 (x) is a sum of four terms,

c 6 (x) = m 6 (x) -Γ (4,2) (x) -Γ (3,3) (x) + 2Γ (2,2,2) (x), (68)
Γ (4,2) (x) = π∈P (4,2) ([[1,6]]) T π . (70)
In the above equation, each T π is computed so that factors of order 4 come first (e.g. there will be no such term as m 2 (x) ⊗ m 4 (x)). For instance,

T {{1;3;5;6};{2;4}} = ( j 1 ,..., j 6 )∈[[1,p]] 6
mom (4) (x j 1 , x j 3 , x j 5 , x j 6 )mom (2) (x j 2 , x j 4 )e j 1 e ′ j 2 ⊗ e j 3 e ′ j 4 ⊗ e j 5 e ′

j 6 = ( j 1 ,..., j 6 )∈[[1,p]] 6 mom (4) (x j 1 , x j 3 , x j 5 , x j 6 )mom (2) (x j 2 , x j 4 )(e j 1 ⊗ e j 3 ⊗ e j 5 )(e ′ j 2 ⊗ e ′ j 4 ⊗ e ′ j 6 ) = ( j 1 ,..., j 6 )∈[[1,p]] 6 mom (4) (x j 1 , x j 3 , x j 5 , x j 6 )mom (2) (x j 2 , x j 4 )(e j 1 ⊗ e j 3 ⊗ e j 5 )(e ′ j 6 ⊗ e ′ j 2 ⊗ e ′ j 4 )P 2,1 P 3,2 = ( j 1 ,..., j 6 )∈[[1,p]] 6
mom (4) (x j 1 , x j 3 , x j 5 , x j 6 )mom (2) (x j 2 , x j 4 )(e j 1 e ′ j 6 ⊗ e j 3 ⊗ e j 5 ) ⊗ (e ′ j 2 ⊗ e ′ j 4 )P 2,1 P 3,2

= M 4 (x) ⊗ vec ′ m 2 (x)P 2,1 P 3,2 . (71)
Following similar lines, T {{1;2;3;4};{5;6}} = m 4 (x) ⊗ m 2 (x) (72) T {{1;2;3;5};{4;6}} = M 4 (x) ⊗ vec ′ m 2 (x) (73) T {{1;2;4;5};{3;6}} = P 3,2 m 4 (x) ⊗ m 2 (x) (74)

T {{1;3;4;5};{2;6}} = M 4 (x) ⊗ vec ′ m 2 (x)P 2,1 (75)
T {{2;3;4;5};{1;6}} = P 2,1 P 3,2 m 4 (x) ⊗ m 2 (x) (76) T {{1;2;3;6};{4;5}} = m 4 (x) ⊗ m 2 (x)P 3,2 (77)
T {{1;2;4;6};{3;5}} = M 4 (x) ′ ⊗ vecm 2 (x) (78) T {{1;3;4;6};{2;5}} = m 4 (x) ⊗ m 2 (x)P 3,2 P 2,1 (79)
T {{2;3;4;6};{1;5}} = P 2,1 M 4 (x) ′ ⊗ vecm 2 (x)
(80) T {{1;2;5;6};{3;4}} = P 3,2 m 4 (x) ⊗ m 2 (x)P 3,2 (81) T {{2;3;5;6};{1;4}} = P 2,1 P 3,2 m 4 (x) ⊗ m 2 (x)P 3,2 (82) T {{1;4;5;6};{2;3}} = P 3,2 m 4 (x) ⊗ m 2 (x)P 3,2 P 2,1 (83) T {{2;4;5;6};{1;3}} = P 3,2 P 2,1 M 4 (x) ′ ⊗ vecm 2 (x) (84) T {{3;4;5;6};{1;2}} = P 2,1 P 3,2 m 4 (x) ⊗ m 2 (x)P 3,2 .P 2,1 .

(

The last step is to collect terms in m 4 (x) ⊗ m 2 (x), M 4 (x) ⊗ vec ′ m 2 (x), and M 4 (x) ′ ⊗ vecm 2 (x) respectively:

Γ (4,2) (x) = (I p 3 + P 3,2 + P 2,1 P 3,2 ) R m 4 (x) ⊗ m 2 (x) + (I p 3 + P 3,2 + P 2,1 P 3,2 ) R m 4 (x) ⊗ m 2 (x)P 3,2 P 2,1 + (I p 3 + P 3,2 + P 2,1 P 3,2 ) R m 4 (x) ⊗ m 2 (x)P 3,2 +M 4 (x) ⊗ vec ′ m 2 (x) (I p 3 + P 2,1 + P 2,1 P 3,2 ) P ′ + (I p 3 + P 2,1 + P 3,2 P 2,1 ) P M 4 (x) ′ ⊗ vecm 2 (x) = R m 4 (x) ⊗ m 2 (x) R ′ + M 4 (x) ⊗ vec ′ m 2 (x) P ′ + P M 4 (x) ′ ⊗ vecm 2 (x) . ( 86
)
Expression of Γ (3,3) (x):

It is seen that P (3,3) ([ [START_REF] Hald | The early history of the cumulants and the gram-charlier series[END_REF][START_REF] Albera | Blind identification of over-complete mixtures of sources (biome)[END_REF]]) contains 6! 2!3!3! = 10 elements:

P (3,3) ([[
Γ (3,3) (x) = π∈P (3,3) ([[1,6]]) T π . (88)
In the above equation, each T π is computed so that factor m 3 (x) comes first (e.g. there will be no such term as m 3 (x) ′ ⊗ m 3 (x)):

T {{1;2;3};{4;5;6}} = m 3 (x) ⊗ m 3 (x) ′ (89) T {{1;2;4};{3;5;6}} = P 2,1 P 3,2 m 3 (x) ⊗ m 3 (x) ′ P 2,1 P 3,2 (90) T {{1;3;4};{2;5;6}} = m 3 (x) ⊗ m 3 (x) ′ P 2,1 (91) T {{2;3;4};{1;5;6}} = P 3,2 m 3 (x) ⊗ m 3 (x) ′ P 2,1 P 3,2 (92) T {{1;2;5};{3;4;6}} = P 3,2 m 3 (x) ⊗ m 3 (x) ′ (93) T {{1;3;5};{2;4;6}} = vecm 3 (x) ⊗ vec ′ m 3 (x) (94) T {{2;3;5};{1;4;6}} = P 2,1 P 3,2 m 3 (x) ⊗ m 3 (x) ′ (95) T {{1;4;5};{2;3;6}} = P 3,2 m 3 (x) ⊗ m 3 (x) ′ P 2,1 (96) T {{2;4;5};{1;3;6}} = m 3 (x) ⊗ m 3 (x) ′ P 2,1 P 3,2 (97) T {{3;4;5};{1;2;6}} = P 2,1 P 3,2 m 3 (x) ⊗ m 3 (x) ′ P 2,1 .

(98) Therefore,

Γ (3,3) (x) = (I p 3 + P 3,2 + P 2,1 P 3,2 ) R m 3 (x) ⊗ m 3 (x) ′ + (I p 3 + P 3,2 + P 2,1 P 3,2 ) R m 3 (x) ⊗ m 3 (x) ′ P 2,1 + (I p 3 + P 3,2 + P 2,1 P 3,2 ) R m 3 (x) ⊗ m 3 (x) ′ P 2,1 P 3,2 + vecm 3 (x) ⊗ vec ′ m 3 = R m 3 (x) ⊗ m 3 (x) ′ P ′ + vecm 3 (x) ⊗ vec ′ m 3 . (99
Γ (2,2,2) (x) = π∈P (2,2,2) ([[1,6]]) T π . (101)
In the above equation, each T π is computed by imposing the ordering m 2 (x) < vecm 2 (x) < vec ′ m 2 (x), which amounts to placing first blocks with polarities {even}{odd}, then {odd}{odd}, and finally {even}{even}:

T {{1;2};{3;4};{5;6}} = m 2 (x) ⊗ m 2 (x) ⊗ m 2 (x) (102)
= P 2,1 m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x) (107) T {{1;2};{4;5};{3;6}} = P 3,2 m 2 (x) ⊗ m 2 (x) ⊗ m 2 (x) (108) T {{1;4};{2;5};{3;6}} = P 3,2 P 2,1 m 2 (x) ⊗ m 2 (x) ⊗ m 2 (x) (109) T {{2;4};{1;5};{3;6}} = P 2,1 m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x)P 2,1 P 3,2 (110) T {{1;3};{4;5};{2;6}} = P 3,2 P 2,1 m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x)P 2,1 (111) T {{1;4};{3;5};{2;6}} = m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x)P 2,1 (112) T {{3;4};{1;5};{2;6}} = P 2,1 m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x)P 2,1 (113) T {{2;3};{4;5};{1;6}} = P 2,1 P 3,2 m 2 (x) ⊗ m 2 (x) ⊗ m 2 (x) (114) T {{2;4};{3;5};{1;6}} = m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x)P 2,1 P 3,2 (115) T {{3;4};{2;5};{1;6}} = P 2,1 P 3,2 P 2,1 m 2 (x) ⊗ m 2 (x) ⊗ m 2 (x). ( 116
P m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x)P 2,1 P 3,2 = Q m 2 (x) ⊗ m 2 (x) ⊗ m 2 (x) + P m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x) P ′ . (117)
This completes the proof of theorem 2.

Two immediate consequences

4.1. Upper bounds on the rank of cumulant matrices Corollary 3. Let x be a real-valued random vector of dimension p × 1 such that E [|x| n ] < ∞ for a given integer n ≥ 6. Then, its cumulant matrices of orders 2 to 6 have their rank upper-bounded as follows:

R c 2 (x) ≤ p (118) R c 3 (x) ≤ p (119) R c 4 (x) ≤ p(p + 1) 2 (120) R c 5 (x) ≤ p(p + 1) 2 (121) R c 6 (x) ≤ p(p + 1)(p + 2) 6 . (122)
Proof. The first two inequalities are obvious. For the other ones, let us introduce

K = 1 2 K and Q = 1 6 Q, (123)
where K and Q are matrices of dimension p 2 × p 2 and p 3 × p 3 , as given in Eqs. (32) and (35), respectively. By using the symmetry properties of cumulants, it comes that

c 4 (x) = Kc 4 (x), c 5 (x) = c 5 (x) K, et c 6 (x) = Qc 6 (x), (124)
from which the following majorations immediately follows:

R c 4 (x) ≤ R K , R c 5 (x) ≤ R K , and R c 6 (x) ≤ R Q . ( 125
)
Since matrices K and Q are idempotent (see Schott [START_REF] Schott | Kronecker product permutation matrices and their application to moment matrices of the normal distribution[END_REF]Theorem 1]), their rank is returned by their trace. Thus,

R K = Tr( K) = 1 2 Tr(I p 2 + K pp ) = 1 2 p 2 + Tr(K pp ) . (126)
The trace of the commutation matrices is in given Magnus and Neudecker [START_REF] Magnus | The commutation matrix: Some properties and applications[END_REF] as

Tr(K ab ) = 1 + gcd(a -1, b -1); (127)
where gcd stands for "greatest common divisor". This implies, in particular, Tr(K pp ) = p; therefore,

R K = 1 2 p 2 + p = p(p + 1) 2 . (128)
Similarly,

R Q = Tr( Q) = 1 6 Tr(Q) = 1 6 Tr I p 3 + I p ⊗ K pp + K pp ⊗ I p + K pp 2 + K p 2 p + (K pp ⊗ I p )K p 2 p = 1 6 p 3 + p 2 + p 2 + p + p + p 2 = p 3 + 3p 2 + 2p 6 = p(p + 1)(p + 2) 6 . (129)

The sixth-order moment matrix of a Gaussian vector

The expression of the sixth-order moment matrix of a Gaussian vector has been the aim of several research works. It has been provided by Neudecker [START_REF] Neudecker | Some theorems on matrix differentiation with special reference to kronecker matrix products[END_REF], Magnus and Neudecker [START_REF] Magnus | The commutation matrix: Some properties and applications[END_REF], Tracy and Sultan [START_REF] Tracy | Higher order moments of multivariate normal distribution using matrix derivatives[END_REF], and more recently by Schott [START_REF] Schott | Kronecker product permutation matrices and their application to moment matrices of the normal distribution[END_REF]. It happens here to be a direct corollary of theorem 2. where matrices Q and P are defined in theorem 2.

Proof. Since z is Gaussian, all its cumulants of order greater than two are nil. Therefore, the only non-zero terms in Eq. ( 12), with k = 6, corresponds to the partition of type (2, 2, 2). Hence,

m 6 (z) = Γ (2,2,2) (x), (131)
where

Γ (2,
This quantity is given by Eq. ( 117), wherein m 2 (x) = c 2 (x) = Σ.

Conclusion

This paper has introduced a methodology for deriving explicit formulae for the cumulant matrices of a random vector in terms of moment matrices of lower order. The methodology reproduces known formulae, up to order four, and returns new ones at order five and six, which are all embodied in theorem 2. In principle, it can also be applied to address higher orders. Such formulae of the cumulant matrices have several advantages. First, from the notational point of view, they are quite concise. Second, from the theoretical point of view, they evidence a direct link with moment matrices. As a consequence, the moment matrices of a random Gaussian vector -which have been of some interest in the literature -are easily found be setting the cumulant matrix to zero (e.g. see corollary ( 4)). In addition, they make possible to upper-bound the rank of the cumulant matrices, as described in corollary (3)). These properties have considerable importance in applications as such independent component analysis (ICA) [START_REF] Comon | Independent component analysis, a new concept?[END_REF], which strongly rely on the use of the cumulant matrices of order four -the Quadricovariance [START_REF] Albera | Icar: a tool for blind source separation using fourth-order statistics only[END_REF] -or six -the Hexacovariance Albera et al. [START_REF] Albera | Blind identification of over-complete mixtures of sources (biome)[END_REF]. In this case, the upper-bound on the rank of the cumulant matrices indicates the maximum number of independent sources that can be separated. Finally, the proposed formulae have a definite advantage when they come to be coded with high-level matrix/array programming language, such as Matlab R . By way of an example, it has been verified by the authors that the numerical computation of the cumulant matrices of orders five and six with dimension p = 6 are about 150 times faster than with coded with Leonov and Shiryaev formula (version 7 of Matlab used on a laptop computer with 2.10GHz clock and 4Go RAM).

  with P s ([[1, 2k]]) the set of partitions of [[1, 2k]] of type s. For instance, for k = 2,P(4) = {(4); (3, 1); (2, 2); (2, 1, 1); (1, 1, 1, 1)} ,(43)and for s = (2, 2), |s| = 2 and P s ([[START_REF] Hald | The early history of the cumulants and the gram-charlier series[END_REF][START_REF] Albera | Icar: a tool for blind source separation using fourth-order statistics only[END_REF]]) = {{{1; 2}; {3; 4}}; {{1; 3}; {2; 4}}; {{1; 4}; {2; 3}}} .

j 4 T

 4 {{1;4};{2;3}} .

  ) T {{1;3;5};{2;4}} = vecm 3 (x)vec ′ m 2 (x) (63) T {{1;4;5};{2;3}} = P 3,2 m 3 (x) ⊗ m 2 (x)K pp (64) T {{2;4;5};{1;3}} = P 3,2 P 2,1 m 3 (x) ′ ⊗ vecm 2 (x) (65) T {{3;4;5};{1;2}} = P 2,1 P 3,2 m 3 (x) ⊗ m 2 (x).K pp .

  T {{1;3};{2;4};{5;6}} = P 3,2 P 2,1 m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x)P 2,1 P 3,2 (103) T {{2;3};{1;4};{5;6}} = P 2,1 m 2 (x) ⊗ m 2 (x) ⊗ m 2 (x) (104) T {{1;2};{3;5};{4;6}} = m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x) (105) T {{1;3};{2;5};{4;6}} = P 3,2 P 2,1 m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x) (106) T {{2;3};{1;5};{4;6}}

Corollary 4 .

 4 Let z ∼ N(0 p , Σ) be a zero-mean, real-valued, Gaussian vector. Its sixth-order moment matrix readsm 6 (z) = E(zz ′ ⊗ zz ′ ⊗ zz ′ ) = Q (Σ ⊗ Σ ⊗ Σ) + P Σ ⊗ vecΣ ⊗ vec ′ Σ P ′(130)

( j 1 2 ⊗ e j 3 e ′ j 4 ⊗ e j 5 e ′ j 6 

 1246 ,..., j 6 )∈[[1,p]] 6 B∈π cum (2) (x j l : l ∈ B)e j 1 e ′ j         .

  = (I p 3 + P 2,1 + P 3,2 + P 3,2 P 2,1 + P 2,1 P 3,2 + P 2,1 P 3,2 P 2,1 )+ (I p 3 + P 2,1 + P 3,2 P 2,1 ) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x)+ (I p 3 + P 2,1 + P 3,2 P 2,1 )

	m 2 (x) ⊗ vecm 2 (x) ⊗ vec ′ m 2 (x)P 2,1
	P
	+ (I p 3 + P 2,1 + P 3,2 P 2,1 )
		)
	Therefore,
	Γ (2,2,2) (x) Q	m 2 (x) ⊗ m 2 (x) ⊗ m 2 (x)

P m 2 (x)


Cuong Kieu
email: vkieu@labri.fr

Muriel Visani
email: muriel.visani@univ-lr.fr

Nicholas Journet
email: journet@labri.fr

Rémy Mullot
email: remy.mullot@univ-lr.fr

Jean-Philippe Philippe Domenger
email: domenger@labri.fr

V C Kieu

An Efficient Parametrization of Character Degradation Model for Semi-synthetic Image Generation

Keywords: degradation model, synthetic document image, degradation model validation, performance evaluation

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

In recent years, degradation models are widely used to generate a benchmarking database in order to assess the performance of different document analysis and recognition methods (e.g. symbol recognition & spotting system evaluation in [START_REF] Delalandre | Generation of Synthetic Documents for Performance Evaluation of Symbol Recognition & Spotting Systems[END_REF], handwriting recognition algorithm in [START_REF] Mori | Generating New Samples from Handwritten Numerals Based on Point Correspondence[END_REF], OCR, text line detection). The other interest of using degradation model for generating synthetic (or semi-synthetic) images is to enrich training databases. For example, the authors Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from Permissions@acm.org. HIP '13 August 24 2013, Washington, DC, USA. Copyright is held by the owner/author(s). Publication rights licensed to ACM. ACM 978-1-4503-2115-0/13/08 ...$15.00. http://dx.doi.org/10.1145/2501115.2501127 of [START_REF] Varga | Effects of Training Set Expansion in Handwriting Recognition Using Synthetic Data[END_REF] and [START_REF] Varga | Generation of Synthetic Training Data for an HMM-based Handwriting Recognition System[END_REF] have improved their system performance by creating synthetic handwritten characters.

Synthetic document image generation is a topic that has been little discussed in the state of the art. Baird presented a state of the art of several methods in [START_REF] Baird | The State of the Art of Document Image Degradation Modeling[END_REF]. He mentioned several degradations that appear in real documents. He also proposed a general methodology for creating degradation models. In [START_REF] Loce | Halftone Banding due to Vibrations in A Xerographic Image Bar Printer[END_REF], the authors proposed a model for reproducing noise that probably appears during the barcode printing process. J. Zhai et al. [START_REF] Jian Zhai | A Line Drawings Degradation Model for Performance Characterization[END_REF] have proposed a model named "hard pencil noise". This model can synthesize white specks to evaluate their line detection method. A modeling of the bleed-through defect has been proposed by R. F. Moghaddam et al. [START_REF] Moghaddam | Low Quality Document Image Modeling and Enhancement[END_REF]. It can reproduce the apparition of recto ink on verso document side background. In [START_REF] Smith | Modeling Image Degradations for Improving OCR[END_REF], the author proposed a noise model that reproduces noise appearing when scanner sensor fails. At last, Kanungo model presented in [START_REF] Kanungo | Global and Local Document Degradation Models[END_REF] (and validated in [START_REF] Kanungo | A statistical, Nonparametric Methodology for Document Degradation Model Validation[END_REF], [START_REF] Li | Validation of Image Defect Models for Optical Character Recognition[END_REF]) can generate the salt and pepper noise that appears frequently near characters. As mentioned by [START_REF] Li | Validation of Image Defect Models for Optical Character Recognition[END_REF] , the difficulty in creating a degradation model is to validate it (i.e. to prove that a model generates realistic defects). This article presents a new contribution started in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF] and [START_REF] Kieu | Semi-synthetic Document Image Generation Using Texture Mapping on Scanned 3D Document Shapes[END_REF] where we proposed two new degradation models. These models are able to reproduce specific defects appearing in old document images. In [START_REF] Kieu | Semi-synthetic Document Image Generation Using Texture Mapping on Scanned 3D Document Shapes[END_REF], we presented a 3D model for reproducing geometric distortions such as folds, torns, or convexo-concaves of the paper sheet. In [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF], we investigated a model which allows to create gray-level defects such as dark specks near characters or ink discontinuities. We are then collaborating with other researchers in order to estimate initially the impact of defects. Recently, 6000 synthetic degraded images for the music score removal staff line competition of ICDAR 20131 were generated by using the two models.

As mentioned earlier in the state of the art, the parametrization of a degradation model is a complex task. First, in order to use a model, it is necessary to understand how to set the different parameters resulting in a lot of possible combina-tions. Then, for specific purpose, an user might have to generate specific defects. For example testing the robustness of an OCR might require to generate characters with many disconnections, whereas testing a word spotting algorithm requires the generation of ink spots between words or lines. Hence, we propose a parametrization method for estimating parameters of the model presented in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF] according to the quantity and the specificity of character defects given by users. Real ancient document images suffer from defects that mostly appear in the neighborhood of the characters. Usually, these noise regions can be divided into three types: independent spot, overlapping spot, and disconnection spot. The independent spots are the noise regions that appear inside or outside of characters. An independent black spot is presented in Fig. 1-a while a white one can be seen in Fig. 1b). The overlapping spots are the noise regions that overlap a character and therefore modify the edge of the character. For example, Fig. 1-c-d have respectively overlapping black spot and overlapping white spot. At last, the disconnection spots are only the white noise regions that break the connectivity of characters (i.e. the disconnection white spot in Fig. 1-e). This paper is organized as follows: first, in Section 2, we explain shortly the principle of our degradation model presented in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF]. Then, in Section 3, an overview of the proposed method for parametrizing automatically our model is described. Several algorithms (for selecting noise regions and type of defect) are introduced in section 4. Finally, experimental results and conclusion are given in section 5 and 6, respectively.

PRESENTATION OF THE CHARACTER DEGRADATION MODEL

As detailed in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF], we proposed a character degradation model based on three steps. Firstly, the centres of noise regions are defined as the flipped pixels resulting from the non-linear local selection process presented by Kanungo et al. [START_REF] Kanungo | Global and Local Document Degradation Models[END_REF]. This process flips, on a binary image, the value of some foreground and background pixels. These pixels, called "seedpoints", are mostly near the edge of characters. They are divided into two sets. The first set (P f b ) represents the centre of each future white spot. The second one (P bf ) represents the centre of each future black spot. This degradation model has five parameters, θ = (α0, α, β0, β, η(η0, η1)) for controlling the amount of generated seed-points.

Secondly, for each seed-point, three properties are computed: shape, direction, and size. In order to create synthetic images as realistic as possible, we propose to generate an elliptic noise region (defined by a minor and a major axis).

The major axis depends on the gradient vector at the centre of this noise region and a parameter a0. The parameter a0 indicates that the generated noise region becomes large or not. The minor axis is set according to the major axis value and a parameter g (i.e. flattening factor of an ellipse). For each noise region, according to the value of the parameters a0 and g, it will be an independent spot, an overlapping spot or a disconnection spot. For example, Fig. 2 provides six degraded images with the three types of degradations. Finally, the gray value of pixels inside each elliptic noise region is changed. To produce a noise as realistic as possible, a new grayscale value is randomly set by using a generation function that satisfies the normal distribution. The mean of this distribution is set according to the distance between each pixel and the centre of the noise region. The variance σ of the function is an input parameter. Thus, our grayscale character degradation model lies on eight parameters: ∆(α0, α, β0, β, η(η0, η1), a0, g, σ) that might lead to difficulties in generating specific synthetic images.

To simplify the noise generation process and to be able to generate specific character defects, we propose a solution for users to parametrize the model with only 4 values: the desired amount of noise regions and the proportion of wished independent spot, overlapping spot, and disconnection spot in the final degraded image.

PROPOSED METHOD OVERVIEW

The noise model presented in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF] can generate various kinds of differently visual artifacts (white/dark spots, character disconnections . . . ). However, with this model, the amount of each type of noise is randomly generated, which makes the visual appearance of the semi-synthetic image result quite unpredictable, and unsuitable to the user's requirement. In this section, we therefore introduce an overview of the proposed method which allows the user to select the amount of each kind of noise injected into the degraded image.

By only using four parameters, the three steps of the degradation model presented in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF] are now simplified as in Fig. 3. First, according to the desired amount of noise regions, we propose an algorithm for setting automatically the Kanungo model parameters [START_REF] Kanungo | Global and Local Document Degradation Models[END_REF] and generating the Nsp seed-points.

Then, with the given number of independent spot, overlapping spot and disconnection spot, we propose a heuristic rule that assigns each seed-point to one of the three categories. Finally, the size of each noise region is set according to another heuristic rule based on a randomly pseudo process. i be the sizes of the major axis of these spots. Fig. 4-a shows that the white spot (with its centre inside a character) is connected with the nearest edge of the stroke when its size ai is equal to a01i. It touches the second one when its size ai is equal to a02i. Indeed, this white spot is an independent white spot until ai = a01i. It becomes an overlapping white spot with a01i ≤ ai ≤ a02i whereas it will be a disconnection spot with ai > a02i. In the same way, the black spot (its centre near a character) C bi will be an independent black spot if a ′ i < a ′ 01i and an overlapping black spot if a ′ 01i ≤ a ′ i ≤ a ′ 02i (see Fig. 4-b) . As mentioned earlier in Section 1, we consider that only white spots can generate disconnection spots. Finally, each seed-point has two thresholds a01i and a02i which will be used to specify the type of each seed-point.

Based on this analysis, the noise classification can be generalized by Definition 3.1.

Definition 3.1 For each seed-point Ci, let a01i, a02i be its thresholds and ai be its major axis size:

• If ai < a01i, Ci is an independent spot.

• If a01i ≤ ai ≤ a02i, Ci is an overlapping spot.

• If ai > a02i and Ci ∈ P f b , Ci is a disconnection spot.

PROPOSED METHOD 4.1 Seed-points position selection

The seed-point generation process detailed in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF] is shortly described in the following procedure (where n is the number of obtained seed-points):

Algorithm 4.1: Seed-PointsGeneration(α 0 , α, β 0 , β, η 0 , η 1 )

n ← 0 for each pixel ∈ Binarized Image do              d ← distance(pixel to character boundary) if (pixel is Background) then weight ← (β 0 e -βd + η 0 ) else weight ← (α 0 e -αd + η 1 ) r ← unif orm random(0, 1) if (weight ≥ r) then pixel ← seed point; n ← n + 1; return (n)
Instead of using the probability procedure, which takes into account the five parameters of Kanungo noise model θ = (α0, α, β0, β, η(η0, η1)), an adaptive procedure is proposed. For this procedure, a user only needs to indicate the desired number of wished seed-points (Nsp) and the wished percentage of disconnection spots (D) among the Nsp seedpoints. As mentioned earlier in [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF], if we generate too many seed-points very close to each other, the result might look too much synthetic. For that reason, we propose to set the maximum number of seed-points that a user can generate to the number of connected components. This value is usually close to the number of characters appearing in the document image. The proposed procedure is detailed as follows: The convergence of Algorithm 4.2 is proved in Appendix A. For the first iteration, α and β are set to 0. It means that the probability for each pixel to be a seed-point is maximum (see ( 3) and (4) in Appendix A). Then, to optimize the while loop, α and β are selected by using a random number generator whose mean is calculated as in (9) (i.e. the expectation E in ( 9) is equal to the number of input seed-points). The parameter η0 is equal to (-D) in order to ensure that we always generate D percent of disconnection spots. For example, if D is set at 1 (100%) by the user, according to (3) and ( 4), the probability for a background pixel to become an ink pixel is equal to (β0e -βd -1) < 0. It means that only disconnection spots will be generated.

Noise region classification

The aim of this step is to assign each seed-point to one of the three types. The difficulty is to find an allocation which generates a realistic semi-synthesized document image. Indeed, if too many large noise regions are generated, it might lead to non-realistic results. For example, if a seed-point far from a character boundary is assigned as an overlapping spot, it will produce a non-realistic big black spot. That is the reason why we propose a heuristic rule which assigns, for each seed-point, a type of noise based on its two thresholds.

The output of the previous step is a set Nsp of seed-points. Let I, O and D be the number of wished independent spots, overlapping spots, and disconnection spots (I + O + D = Nsp). Let Nis, Nos, N ds be the sets of independent spots, overlapping spots, and disconnection spots. Let min_a01 (setOfSeedPoint S) and min_a02 (setOfSeedPoint S) be the two functions which respectively return the seed-point with the lowest a01 and a02 value. Let's see an example of five seed-points (C1 → C5) given in Table 1. The percentages of the three types are given by a user as follows: 20% of independent spots, 40% of overlapping spots and 40% of disconnection spots. Indeed, with I = 1, O = 2 and D = 2, Nis will have one seed-point, Nos and N ds will have two seed-points after the application of Algorithm 4.1. The two values (a01 and a02 ) are computed for each of the five seed-points. First, C1 and C4 are set as disconnection spot. Then, C2 and C3 are set as Overlaping spot. At last, C5 is set as independent spot.

Noise generation process

As a result of the previous step, all the seed-points are classified according to the user's wishes. In this step, the value of every pixel inside each ellipse is degraded to generate noise regions. In order to generate degraded images as realistic as possible, a heuristic rule is used to set randomly the size of the major axis of each ellipse according to the type of degradation while the Definition 3.1 is respected. Let ai, aj, and a k be respectively the sizes of the major axis of an elliptic independent spot, overlapping spot, and disconnection spot regions. Therefore:

   ai = a01i × µi, 0 < µi < 1 aj = a01j + µj × (a02j -a01j), 0 < µj < 1 a k = a 02k + µ k × δ, 0 < µ k < 1 (1)
To define these sizes, the values of µi, µj, µ k are automatically selected by using a random number generator ranging from 0 to 1. The flattening factor g of the elliptic noise region can be randomly chosen between [0, 1] as well. The parameter δ is set by the average width of all the connected components to avoid generating large disconnection spots. The pixel values within the elliptic noise region are modified by a random number generator following the normal distribution as presented in Section 2.3 of [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF]. Last, the Gaussian blur operation is appliqued on the noise region so as to make the values of two adjacent pixels more coherent.

Degradation level estimation

As a method using a random process for generating degraded document images, a value that indicates the effected proportion of degradation of the image result needs to be computed. In general, this value, also known as the degradation level, is calculated by comparing the image before and after the random degradation process. In [START_REF] Jian Zhai | A Line Drawings Degradation Model for Performance Characterization[END_REF], the authors calculated the value as the difference of gray value of all pixels before and after the degradation process. Let L be the difference between the image before and after this process. Let I and I ′ be the gray images before and after the process, respectively. Therefore, L can be calculated as follows:

L = I -I ′ (2)

EXPERIMENTAL RESULTS

Generation of images with a fixed number of seed-points

In this test, the number of seed-points is fixed by the number of connected components. Fig. 5 provides three degraded images of the original one in Fig. 5-a, which has 351 connected components. The image in Fig. 5-b is degraded by generating only independent spots. It has the effect of generating only small black noise regions near the characters or small white ones inside the characters. In Fig. 5-c, with the same number of seed-points, the degraded image contains only overlapping spots in which black overlapping spots connect with the edge of characters whereas white ones modify the ink character regions. In Fig. 5-d, we decided to generate only disconnection spots. Indeed, a lot of strokes of characters are completely disconnected. Moreover, the degradation values (L) of the three images are coherent with a visual validation. For the three images, the degradation values are L = 183, L = 1623 and L = 4026, respectively. This shows that the document image in Fig. 5-d is the most degraded image of the three images. Fig. 5.e h present three zoomed degraded versions of the word "que" at the end of the text line 9 of the three images.

Generation of images with different amount of seed-points

In this second test, four degraded images are generated from the original image presented in Fig. 6 with the fixed proportion of the three types of degradations (I = 15%, O = 60%, D = 25%) whereas the number of wished seed-points increases progressively.

Fig. 7 provides the degraded letters of the four images which are visually realistic in spite of the high density of seedpoints around the characters (e.g. in Fig. 7-e, there are nine

CONCLUSION

In this paper, we present an efficient parametrization method for a character degradation model. This method aims at simplifying the eight model parameters ∆(α0, α, β0, β, η, a0, g, σ) by using three steps. First, Nsp seed-points are selected by adjusting automatically the parameters of Kanungo noise model. The two thresholds a01, a02 of each seed-point are also estimated at this step. Then, the seed-points can classified into three types of degradation. Finally, the size of each noise region is randomly chosen in order to generate a specific noise region. Therefore, the model only takes into account four parameters: the number of seed-points and the percentages of the three types of degradation that a user wants to generate. Many images have been generated by using our method. All the images satisfy the user's wishes and look realistic. Since, the calculated degradation levels are coherent with the visual results, that is encouraging us to estimate automatically the four parameters by only the degradation level given by user. The formally validation of the model is also our future work.

Let Igray be the input image. Let n be the number of all pixels in Igray and Z be the maximum distance of a pixel to the nearest edge of characters (the distance transform). The weight w of each pixel P (x,y) is calculated as follows:

w =   
α0 × e -αd + η1, P : F oreground β0 × e -βd + η0, P : Background

(d ∈ [1 → Z]) (3)
where d is the distance transform of the pixel P . As mentioned earlier in Algorithm 4.1, a pixel will be a seed-point when its weight is greater than or equals to the random number r. The random number r is generated by an uniform distribution number generator. Let's see a pixel P (x,y) which has the weight wi. The probability of P (x,y) is calculated hereafter (see Fig. 8) so that this pixel will be a seed-point: According to (4), the probability with which k pixels will be seed-points is calculated as follows:

pi = wi (4)
P (nsp=k) = m=M Im⊂[1...n];|Im|=k;l=1 ( i / ∈Im (1 -wi) × j∈Im wj) (5)
where M is the number of subsets Im ⊂ [1 . . . n] in which each subset has k elements. Thus, M is calculated as follows:

M = n k = n! k!(n -k)! (6)
According to (3), ∀ pixel P: 0 < wi < 1; therefore, 0 < P (nsp=k) < 1. This means that we can always generate k seed-points with the probability calculated as in [START_REF] Kanungo | Global and Local Document Degradation Models[END_REF]. As a result, the Algorithm 4.2 always converges. In addition, the expected value is calculated as follows:

E = n l=0 l × p (k=l) = n i=1 wi (7)
Equation [START_REF] Kieu | A Character Degradation Model for Grayscale Ancient Document Images[END_REF] shows that if α and β increase, the weight of each pixel will decrease exponentially according to (3); therefore, the expected number of seed-points will also decrease exponentially.

Let ni1 be the number of Foreground pixels which have the same distance di. Let ni0 be the number of Background pixels which have the same distance di. According to (3) and ( 7), the expectation will be:

E = Z i=1 (ni1e -αd i + ni0e -βd i + ni1η1 + ni0η0) = e α Z i=1 ni1e -d i + e β Z i=1 ni0e -d i + Z i=1 (ni1η1 + ni0η0) = e α A + e β B + C (8)
where A = Z i=1 ni1e -d i , B = Z i=1 ni0e -d i , C = Z i=1 (ni1η1+ ni0η0). Note that A, B, C are the constants for each input image. We suppose that α and β are generated by using a random number generator. The mean µ αβ of this generator can be calculated by using ( 8):

Figure 1 :

 1 Figure 1: Examples of three types of noise in real ancient documents: (a), (b) two independent black and white spots; (b), (c) two overlapping black and white spots; (e) disconnection white spot

7 Figure 2 :

 72 Figure 2: Examples of different noise generation with different a0 values: (a), (e) original images; (b), (f ) two independent white and black spots; (c) overlapping white spot and (g), (h) two overlapping black spots; (d) disconnection spot.

Figure 3 :

 3 Figure 3: The three main steps of our automatic parameter process for character degradation model Before describing the three main steps of the model in detail, we would like to focus on what happens to an elliptic noise region when its size increases. Supposed that two seedpoints Cwi ∈ P f b and C bi ∈ P bf are respectively the centres of future white and black spots. Let ai and a ′i be the sizes of the major axis of these spots. Fig.4-a shows that the white spot (with its centre inside a character) is connected with the nearest edge of the stroke when its size ai is equal to a01i. It touches the second one when its size ai is equal to a02i. Indeed, this white spot is an independent white spot until ai = a01i. It becomes an overlapping white spot with a01i ≤ ai ≤ a02i whereas it will be a disconnection spot with ai > a02i. In the same way, the black spot (its centre near a character) C bi will be an independent black spot if a ′ i < a ′ 01i and an overlapping black spot if a ′ 01i ≤ a ′ i ≤ a ′ 02i (see Fig.4-b) . As mentioned earlier in Section 1, we consider that only white spots can generate disconnection spots.

Figure 4 :

 4 Figure 4: The evolution of a noise region when its major axis size ai increases. (a) case of white spot inside a character. (b) case of black spot near a character.

Algorithm 4 . 2 :

 42 Seed-PointsIdentification(nsp, D) α 0 ← β 0 ← 1, α ← β ← 0 η 0 ← (-D), η 1 ← 0 n ← Seed-PointsGeneration(α 0 , α, β 0 , β, η 0 , η 1 ) while (n = nsp) do    α = random[0, +∞] β = random[0, +∞]n ← Seed-PointsGeneration(α 0 , α, β 0 , β, η 0 , η 1 ) return (α, β, α 0 , β 0 , η 0 , η 1 , n)

Algorithm 4 . 3 :

 43 Seed-PointsAffecation(Nsp, I, O, D) N is , Nos, N ds ← {} sp1, sp2 : seedP oints; while (Nsp is not empty and D > 0) do    sp1 ← min a02(Nsp); if (D > 0 and sp1 ∈ P f b ) then D --; Nsp ← Nsp -{sp1} ; N ds ← N ds + {sp1} ; while (Nsp is not empty and O > 0) do    sp2 ← min a01(Nsp); if (O > 0) then O --; Nsp ← Nsp -{sp2} ; Nos ← Nos + {sp2} ; M ove all the rest seedpoints of Nsp to N is ; return (N is , Nos, N ds )

Figure 6 :Figure 7 :

 67 Figure 6: Original document image with 75 connected components

Figure 8 :

 8 Figure 8: Uniform Distribution

E

  = Ae µ αβ + Be µ αβ + C ⇔ µ αβ = ln( E-C A+B ) (9) (a) Original image with 351 connected components (b) Nsp = 351, I = 100%, O = 0%, D = 0%, L = 183 (c) Nsp = 351, I = 0%, O = 100%, D = 0%, L = 1623 (d) Nsp = 351, I = 0%, O = 0%, D = 100%, L = 4026 (e) Zoomed word original (f) Zoomed word of (b) (g) Zoomed word of (c) (h) Zoomed word of (d)

Figure 5 :

 5 Figure 5: Degraded images with the maximum number of seed-points Nsp = Nccs = 351: (a) original image; (b) degraded image with only the independent spot type and L = 183; (c) degraded image with only the overlapping spot type and L = 1623; (d) degraded image with only the disconnection spot type and L = 4026; (e), (f ), (g), and (h) are four zoomed words of the images (a), (b), (c), (d), respectively.

Table 1 :

 1 Example of noise region classification

	Seed-point a 01 a 02 Final classification
	C 1	1.5	2.4 Disconnection spot
	C 2	1.3	3.7 Overlapping spot
	C 3	2.1	4.6 Overlapping spot
	C 4	1.9	2.7 Disconnection spot
	C 5	2.8	5.4 Independent spot

http://www.cvc.uab.es/cvcmuscima /competition2013/index.htm

http://www.labri.fr/perso/vkieu/content /DegradationModels/cdm model.html

Acknowledgement

This research is done within the DIGIDOC project financed by the ANR (Agence Nationale de la Recherche).