Yoan Miche

Patrick Bas

Amaury Lendasse

Christian Jutten

Olli Simula

Reliable Steganalysis Using a Minimum Set of Samples and Features

This paper proposes to determine a sufficient number of images for reliable classification, and to use feature selection to select most relevant features for achieving reliable steganalysis. First dimensionality issues in the context of classification are outlined and the impact of the different parameters of a steganalysis scheme (the number of samples, the number of features, the steganography method and the embedding rate) is studied. On one hand, it is shown that using Bootstrap simulations, the standard deviation of the classifications results can be very important if too small training sets are used; moreover a minimum of 5000 images is needed in order to perform reliable steganalysis. On the other hand, we show how the feature selection process using the OP-ELM classifier enables both to reduce the dimensionality of the data and to highlight weaknesses and advantages of the six most popular steganographic algorithms.

I. INTRODUCTION

Steganography has been known and used for a very long time, as a way to exchange information in an unnoticeable manner between parties, by embedding it in another, apparently innocuous, document.

Nowadays steganographic techniques are mostly used on digital content. The online newspaper Wired News, reported in one of its articles [START_REF] Mccullagh | Secret Messages Come in[END_REF] on steganography that several steganographic contents have been found on web sites with very large image database such as eBay. Niels Provos [START_REF] Provos | Detecting steganographic content on the internet[END_REF] has somewhat refuted these facts by analyzing and classifying two million images from eBay and one million from USENet network and not finding any steganographic content embedded in these images. This could be due to many reasons, such as very low payloads, making the steganographic images less detectable to steganalysis and hence more secure.

In practice the concept of security for steganography is difficult to define, but Cachin in [START_REF] Cachin | An information-theoretic model for steganography[END_REF] mentions a theoretic way to do so, based on the Kullback-Leibler divergence. A stego process is thus defined as ǫ-secure if the Kullback-Leibler divergence δ between the probability density functions of the cover content p cover and of this very same content embedding a message p stego (that is, stego), is less than ǫ: δ(p cover , p stego ) ≤ ǫ.

(

) 1
The process is called secure if ǫ = 0, and in this case, the steganography is perfect, creating no statistical differences by Y.M., A. L and O.S are from Helsinki University of Technology -Laboratory of Computer and Information Science P.O. Box 5400, FI-02015 HUT, FINLAND ; Y.M., P.B and C.J are from GIPSA-lab Images and Signal Department CNRS, INPG, UJF, Grenoble, France the embedding of the message. Steganalysis would then be impossible.

Fortunately, such high performance for a steganographic algorithm is not achievable when the payload (the embedded information) is important enough; also, several schemes have weaknesses.

One way of mesuring the payload is the embedding rate, defined as follows:

Let S be a steganographic algorithm and M be a cover medium. S, by its design, claims that it can embed at most T Max information bits within M ; T Max is called the capacity of the medium and highly depends on the steganographic (stego) algorithm as well as the cover medium itself. The embedding rate T is then defined as the part of T Max used by the information to embed.

For T i bits to embed in the cover medium, the embedding rate is then T = T i /T Max , usually expressed as percentage. There are other ways to measure the payload and the relationship between the amount of information embedded and the cover medium, such as the number of bits per non zero coefficient. Meanwhile, the embedding rate has the advantage of taking into account the stego algorithm properties and is not directly based on the cover medium properties -since it uses the stego algorithm estimation of the maximum capacity. Hence it has been chosen for this analysis of stego schemes. This paper is focused onto classical feature-based steganalysis. Such steganalysis typically uses a certain amount of images for training a classifier: features are extracted from the images and fed to a binary classifier (usually Support Vector Machines) for training. The output of this classifier is "stego" (modified using a steganographic algorithm) or "cover" (genuine). This classical process is illustrated on Fig. 1 for the part without parenthesis.

The emphasis in this paper is more specifically on the issues related to the increasing number of features, linked to the universal steganalyzers. Indeed, the very first examples of LSB-based steganalysis made use of less than ten features, with an adapted and specific methodology for each stego algorithm. The idea of "universal steganalyzers" then became popular. In 1999, Westfeld proposes a χ 2 -based method, on the LSB of DCT coefficients [START_REF] Westfeld | Attacks on steganographic systems[END_REF]. Five years after, Fridrich in [START_REF] Fridrich | Feature-based steganalysis for jpeg images and its implications for future design of steganographic schemes[END_REF] uses a set of 23 features obtained by normalizations of a much larger set, whilst Farid et al. already proposed in 2002 a set of 72 features [START_REF] Lyu | Detecting hidden messages using higher-order statistics and support vector machines[END_REF]. Some feature sets [START_REF] Agaian | New multilevel dct, feature vectors, and universal blind steganalysis[END_REF] also have variable size depending on the DCT block sizes. Since then, an increasing number of research works are using supervised learning based Fig. 1. Overview of the typical global processing for an analyzed image: features are first extracted from the image and then processed through a classifier to decide whether the image is cover or stego. In the proposed processing is added an extra step aimed at reducing the features number and having an additional interpretability of the steganalysis results, by doing a feature selection. classifiers in very high dimensional spaces. The recent work of Y. Q. Shi et al. [START_REF] Shi | A markov process based approach to effective attacking jpeg steganography[END_REF] is an example of an efficient result but using 324 features based on JPEG blocks differences modeled by Markov processes. This short survey of some feature sets for steganalysis is by no means exhaustive.

All these results do achieve better and better performance in terms of detection rate and enable to detect most stego algorithm for most embedding rates. Meanwhile, there are some side-effects to this growing number of features. It has been shown for example in [START_REF] Franc ¸ois | High-dimensional data analysis: optimal metrics and feature selection[END_REF] that the feature space dimensionality in which the considered classifier is trained, can have a significant impact on its performances: a too small amount of images regarding dimensionality (the number of features) might lead to an improper training of the classifier and thus to results with a possibly high statistical variance.

Comparison of steganalysis methods has been recently proposed by Ker [START_REF] Ker | The ultimate steganalysis benchmark?[END_REF], by focusing on the pdf of one output of the classifier. In this paper is addressed the idea of a practical way of comparing steganographic schemes, which requires a study on multiple parameters that can influence the performance:

1) the number of images used during the training of the classifier: How to determine a sufficient number of images for an efficient and reliable classification (meaning that final results have acceptable variance)? 2) the number of features used: What are the sufficient and most relevant features for the actual classification problem? 3) the steganography method: Is there an important influence of the stego algorithm on the general methodology? 4) the embedding rate used: Does the embedding rate used for the steganography modifies the variance of the results and the retained best features (by feature selection)?

The next section details some of the problems related to the number of features used (dimensionality issues) and commonly encountered in steganalysis: the empty space and the distance concentration phenomena, the important variance of the results obtained by the classifier whenever the number of images used for training is not sufficient regarding the number of features, and finally, the lack of interpretability of the results because of the high number of features. In order to address these issues, the methodology sketched on Fig. 1 is used and more thoroughly detailed: a sufficient number of images regarding the number of features is first established so that the classifier's training is "reliable" in terms of variance of its results; then, using feature selection the interpretability of the results is improved.

The methodology is finally tested in section IV with six differents stego algorithms, each using four different embedding rates. Results are finally interpreted thanks to the most relevant selected features for each stego algorithm. A quantitative study of selected features combinations is then provided.

II. DIMENSIONALITY ISSUES AND METHODOLOGY

The common term "curse of dimensionality" [START_REF] Bellman | Adaptive control processes: a guided tour[END_REF] refers to a wide range of problems related to a high number of features. Some of these dimensionality problems are considered in the following, in relation with the number of images and features.

A. Issues related to the number of images

1) The need for data samples: In order to illustrate this problem in a low-dimensional case, one can consider four samples in a two dimensional space (corresponding to four images out of which two features have been extracted); the underlying structure leading to the distribution of these four samples seems impossible to infer, and so is the creation of a model for it. Any model claiming it can properly explain the distribution of these samples will behave erratically (because it will extrapolate) when a new sample is introduced. On the contrary, with hundreds to thousands of samples it becomes possible to see clusters and relationships between dimensions.

More generally, in order for any tool to be able to analyze and find a structure within the data, the number of needed samples is growing exponentially with the dimensionality. Indeed, consider a d-dimensional unit side hypercube, the number of samples needed to fill the Cartesian grid of step ǫ inside of it is growing as O((1/ǫ) d ). Thus using a common grid of step 1/10 in dimension 10, it requires 10 10 samples to fill the grid.

In practice, most data sets in steganalysis use at least 10 to 20 dimensions, implying a "needed" number of samples impossible to achieve: storing and processing such number of images is currently impossible. As a consequence, the feature space is not filled with enough data samples to estimate the density with reliable accuracy; which can give wrong or high variance models when building classifiers, having to extrapolate for the missing samples: results obtained can have rather high confidence interval and hence be statistically irrelevant. A claim of performance improvement of 2% using a specific classifier/steganalyzer/steganographic scheme with a variance of 2% is rather meaningless.

2) The increasing variance of the results: The construction of a proper and reliable model for steganalysis is also related to the variance of the results it obtains. Only experimental results are provided to support this claim: with a low number of images regarding the number of features (a few hundreds of images for 200 features), the variance of the classifier's results can be very important.

When the number of images increases, this variance decreases toward proper values for classical steganalysis and performances comparisons. These claims are verified in the next section with the experiments.

3) Proposed solution to the lack of images: Overall, these two problems lead to the same conclusion: the number of images has to be important, regarding dimensionality. Theory states that this number is exponential with the number of features, which is impossible to reach for classical steganalysis. Hence, the first step of the proposed methodology is to find a "sufficient" number of images for the number of features used, according to a criterion on the variance of the results.

A Bootstrap [START_REF] Efron | An Introduction to the Bootstrap[END_REF] is proposed for that task: the number of images used for the training of the classifier is increased and for each different number of images, the variance of the results of the classifier is assessed. Once the variance of the classifier is below a certain threshold, a sufficient number of images has been found (regarding the classifier and the feature set used).

B. Issues related to the number of features

1) The empty space phenomenon: This phenomenon was first introduced by Scott and Thompson [START_REF] Scott | Probability density estimation in higher dimensions[END_REF] can be explained with the following example: draw samples from a normal distribution (zero mean and unit variance) in dimension d and consider the probability to have a sample at distance r from the mean of the distribution (zero). It is given by the probability density function:

f (r, d) = r d-1 2 d/2-1 . e -r 2 /2 Γ(d/2) (2)
having its maximum at r = √ d -1. Thus, when dimension increases, samples are getting farther from the mean of the distribution. A direct consequence of this is that for the previously mentionned hypercube in dimension d, the "center" of it will tend to be empty, since samples are getting concentrated in the borders and corners of the cube.

Therefore, whatever model is used in such a feature space will be trained on scattered samples which are not filling the feature space at all. The model will then not be proper for any sample falling in an area of the space where the classifier had no information about during the training. It will have to extrapolate its behavior for these empty areas and will have unstable performances.

2) Lack of interpretability for possible "reverseengineering":

The interpretability (and its applications) is an important motivation for feature selection and dimensionality reduction: high performances can indeed be reached using the whole 193 features set used in this paper for classification.

Meanwhile, if looking for the weaknesses and reasons why these features react vividly to a specific algorithm, it seems rather impossible on this important set.

Reducing the required number of features to a small amount through feature selection enables to understand better why a steganographic model is weak on these particular details, highlighted by the selected features. Such analysis is performed in section IV-C for all six steganographic process. Through the analysis of these selected features, one can consider a "reverse-engineering" of the stego algorithm as illustrated on Fig. 2. By the identification of the most relevant features, the main characteristics of the embedding method can be inferred and the steganographic algorithm be identified if known, or simply understood.

3) Proposed solution to the high number of features: These two issues motivate the feature selection process: if one can reduce the number of features (and hence the dimensionality), the empty space phenomena will have a reduced impact on the classifier used. Also, the set of features obtained by the feature selection process will give insights on the stego scheme and its possible weaknesses.

For this matter, a classical feature selection technique has been used as the second step of the proposed methodology.

The following methodology is different from the one presented previously in [START_REF] Miche | Extracting relevant features of steganographic schemes by feature selection techniques[END_REF], [START_REF] Miche | A feature selection methodology for steganalysis[END_REF]. Indeed, in this article, the goal is set toward statistically reliable results. Also, feature selection has the advantage of reducing the dimensionality of the data (the jumber of features), making the classifier's training much easier. The interpretation of the selected features is also an important advantage (compared to having only the classifier's performance) in that it gives insights on the weaknesses of the stego algorithm.

III. METHODOLOGY FOR BENCHMARKING OF STEGANOGRAPHIC SCHEMES

Addressed problems

The number of data points to be used for building a model and classification is clearly an issue and in the practical case, how many points are needed in order to obtain accurate results -meaning results with small standard deviation.

Reduction of complexity is another main addressed concern in this framework. Then for the selected number of points to be used for classification, and also the initial dimensionality given by the features set, two main steps remain:

• Choosing the feature selection technique. Since analysis and computation can hardly be done on the whole set of features, the technique used to reduce the dimensionality has to be selected.

• Building a classifier; this implies choosing it, selecting its parameters, training and validating the chosen model. The following paragraphs presents the solutions for these two major issues, leading to a methodology combining them, presented on Fig. 3.

A. Presentation of the classifier used: OP-ELM

The Optimally-Pruned Extreme Learning Machine (OP-ELM [START_REF] Miche | A methodology for building regression models using extreme learning machine: OP-ELM[END_REF], [START_REF] Sorjamaa | Long-term prediction of time series using nne-based projection and op-elm[END_REF]) is a classifier based on the original Extreme Learning Machine (ELM) of Huang [START_REF] Huang | Extreme learning machine: Theory and applications[END_REF]. This classifier makes use of single hidden layer feedforward neural networks (SLFN) for which the weights and biases are randomly initialized. The goal of the ELM is to reduce the length of the learning process for the neural network, usually very long (for example if using classical back-propagation algorithms). The two main theorems on which ELM is based will not be discussed here but can be found in [START_REF] Huang | Extreme learning machine: Theory and applications[END_REF]. Fig. 4 illustrates the typical structure of a SLFN (simplified to a few neurons in here). Supposing the neural network is approximating the output Y = (y 1 , . . . , y N ) perfectly, we would have:

M i=1 β i f (w i x j + b i ) = y j , j ∈ 1, N , (3)
with N the number of inputs X = (x 1 , . . . , x N ) (number of images in our case) and M the number of neurons in the hidden layer.

As said, the novelty introduced by the ELM is to initialize the weights W and biases B randomly. OP-ELM, in comparison to ELM, brings a greater robustness to data with possibly dependent/correlated features. Also, the use of other functions f (activation functions of the neural network) makes it possible to use OP-ELM for case where linear components have an important contribution in the classifier's model, for example.

The validation step of this classifier is performed using classical Leave-One-Out cross-validation, much more precise than a k-fold cross-validation and hence not requiring any test step [START_REF] Efron | An Introduction to the Bootstrap[END_REF]. It has been shown on many experiments [START_REF] Miche | A methodology for building regression models using extreme learning machine: OP-ELM[END_REF], [START_REF] Sorjamaa | Long-term prediction of time series using nne-based projection and op-elm[END_REF], that the OP-ELM classifier has results very close to the ones of a Support Vector Machine (SVM) while having computational times much smaller (usually from 10 to 100 times).

B. Determination of a sufficient number of images

A proper number of images, regarding the number of features, has to be determined. Since theoretical values for that number are not reachable, a sufficient number regarding a low enough value of the variance of the results is taken instead (standard deviation will be used instead of variance, in the following).

The OP-ELM classifier is hence used along with a Bootstrap algorithm [START_REF] Efron | An Introduction to the Bootstrap[END_REF] over 100 repetitions; a subset of the complete dataset (10000 images, 193 features) is randomly drawn (with possible repetitions). The classifier is trained with that specific subset. This process is repeated 100 times (100 random drawings of the subset) to obtain a statistically reliable estimation of the standard deviation of the results. The size of the subset drawn from the complete dataset is then increased, and the 100 iterations are repeated for this new subset size.

The criterion to stop this process is a threshold on the value of the standard deviation of the results. Once the standard deviation of the results gets lower than 1%, it is decided that the subset size S giving such low standard deviation of the classifier's results, is sufficient. S is then used for the rest of the experiments as a sufficient number of images regarding the number of features in the feature set.

C. Dimensionality reduction: feature selection by Forward with OP-ELM

Given the sufficient number of images for reliable training of the classifier, S, feature selection can be performed. The second step of the methodology, a Forward algorithm with OP-ELM (Fig. 3) is used.

1) The Forward algorithm: The forward selection algorithm is a greedy algorithm [START_REF] Rossi | Mutual information for the selection of relevant variables in spectrometric nonlinear modelling[END_REF]; it selects one by one the dimensions, trying to find the one that combines best with the already selected ones. The algorithm is (with x i denoting the i-th dimension of the data set)

Algorithm 1 Forward R = {x i , i ∈ 1, d } S = ∅ while R = ∅ do for x j ∈ R do Evaluate performance with S ∪ x j end for Set S = S ∪ {x k }, R = R -x k with x k the

dimension giving the best result in the loop end while

This algorithm requires d(d -1) 2 instances to terminate (to be compared to the 2 d -1 instances for an exhaustivesearch), which might reach the computational limits, depending on the number of dimensions and time to evaluate the performance with one set. With the OP-ELM as a classifier, computational times are not much of an issue. Even if its capacity to isolate efficient features is clear, the forward technique has some drawbacks. First, if two features present good results when put together but poor results if only Bootstrap Fig. 3. Schematic view of the proposed methodology: (1) An appropriate number of data samples to work with is determined using a Bootstrap method for statistical stability; (2) The forward selection is performed using an OP-ELM classifier to find a good features set. Follows a possible interpretation of the features or the typical classification for steganalysis. one of them is selected, forward might not take these into account in the selection process.

Second, it does not allow to "go back" in the process, meaning that if performances are decreasing along the selection process, and that the addition of another feature makes performances increase again, combinations of previously selected features with this last one are not possible anymore.

The Forward selection is probably not the best possible feature selection technique, and recent contribution to these techniques such as Sequential Floating Forward Selection (SFFS) [START_REF] Ververidis | Fast and accurate sequential floating forward feature selection with the bayes classifier applied to speech emotion recognition[END_REF] and improvements of it [START_REF]Fast Sequential Floating Forward Selection applied to Emotional Speech Features estimated on DES and SUSAS Data Collections[END_REF] have shown that the number of computations required for feature selection can be reduced drastically. Nevertheless, the feature selection using Forward has been showing very good results and seems to perform well on the feature set used in this paper. It is not used here in the goal of obtaining the best possible combination of features, but more to reduce the dimensionality and obtain some meaning out of the selected features. Improvements of this methodology could make use of such more efficient techniques of feature selection.

D. General Methodology

To summarize, the general methodology on Fig. 3 uses first a Bootstrap with 100 iterations on varying subsets sizes, to obtain a sufficient subset size and statistically reliable classifiers' results, regarding the number of features used. With this number of images, feature selection is performed using a Forward selection algorithm; this enables to highlight possible weaknesses of the stego algorithm. This methodology has been applied to six popular stego algorithms for testing. Experiments and results, as well as a discussion on the analysis of the selected features are given in the next section.

IV. EXPERIMENTS AND RESULTS

A. Experiments setup 1) Steganographic algorithms used:

Six different steganographic algorithms have been used: F5 [START_REF] Westfeld | F5-a steganographic algorithm[END_REF], Model-Based (MBSteg) [START_REF] Sallee | Model-based steganography[END_REF], MMx [START_REF] Kim | Modified matrix encoding technique for minimal distortion steganography[END_REF] (in these experiments, MM3 has been used), JP Hide and Seek [12], OutGuess [START_REF] Provos | Defending against statistical steganalysis[END_REF], StegHide [START_REF] Hetzl | A graph-theoretic approach to steganography[END_REF]; all of them with four different embedding rates: 5, 10, 15 and 20%.

2) Generation of image database: The image base was constituted of 10000 images from the BOWS2 Challenge [START_REF]Watermarking Virtual Laboratory (Wavila) of the European Network of Excellence ECRYPT[END_REF] database (hosted by Andreas Westfeld [START_REF] Westfeld | Reproducible signal processing (bows2 challenge image database[END_REF]). These images are 512 × 512 PGM greyscale (also available in color).

While the steganographic processes and the proposed methodology for dimensionality reduction and steganalysis are only performed on these 512 × 512 images, any size of image would work just as well.

3) Extraction of the features:

In the end, the whole set of images is separated in two equal parts: one is kept as untouched cover while the other one is stego with the six steganographic algorithms at four different embedding rates: 5%, 10%, 15% and 20%. Fridrich's 193 DCT features [START_REF] Pevny | Merging markov and dct features for multi-class jpeg steganalysis[END_REF] have been used for the steganalysis.

B. Results

Results are presented following the methodology steps. A discussion over the selected features and the possible interpretability of it are developed afterward. In the following, the term "detection rate" stands for the performance of the classifier on a scale from 0 to 100% of classification rate. It is a measure of the performance instead of a measure of error.

1) Determination of sufficient number of samples:

Presented first is the result of the evaluation of a sufficient number of images, as explained in the previous methodology, on Fig. 5. The Bootstrap (100 rounds) is used on randomly taken subsets of 200 up to 9000 images out of the whole 10000 from the BOWS2 challenge.

It can be seen on Fig. 5 that the standard deviation behaves as expected when increasing the number of images, for the cases of JPHS, MBSteg, MMx, OutGuess and StegHide: its value decreases and tends to below 1% of the best performance when the number of images is 5000 (even if for MBSteg with embedding rate of 5% it is a bit above 1%). This sufficient number of samples is kept as the reference and sufficient number. Another important point is that with very low number of images (100 in these cases), the standard deviation is between 1 and about 6.5% of the average classifier's performance; meaning that results computed with as small number of images have at most a ±6.5% confidence interval. While the plots decrease very quickly when increasing the number of images, values of the standard deviation remain very high until 2000 (f) StegHide Fig. 5. Standard deviation in percentage of the average classification result versus the number of images, for all six steganographic algorithms, for the four embedding rates: black circles ( ) for 20%, green squares (2) for 15%, red crosses (×) for 10% and blue triangles ( ) for 5%. These estimations have been performed with the Bootstrap runs (100 iterations). Plots do not have the same scale, vertically.

images; these results have to take into account the embedding rate, which tends to make the standard deviation higher as it decreases.

Indeed, while differences between 15 and 20% embedding rates are not very important on the four previously mentioned stego algorithms, there is a gap between the 5 -10% plots and the 20% ones. This is expected when looking at the performances of the steganalysis process: low embedding rates tend to be harder to detect, leading to a range of possible performances wider than with high embedding rates. Fig. 6 illustrates this idea on the cases of JPHS, MMx, MBSteg StegHide and OutGuess (the F5 results are behaving differently on these plots; this is discussed below).

Finally, the F5 algorithm have "erratic" behaviour, on Fig. 5. This can be explained by the very high performance of the classifier (close to 100%) and hence, the very small possible standard deviation: it ranges between 0.25% and 0.05%. It is then no surprise that the plots have such an erratic behaviour.

The final "sufficient" number of samples retained for the second step of the methodology -the feature selectionis 5000, for two reasons: first, the computational times are acceptable for the following computations (feature selection step with training of classifier for each step); second, the standard deviation is small enough to consider that the final classification results are given with at most 1% of standard deviation (in the case of MBSteg at 5% of embedding rate).

2) Forward feature selection: Features have first been ranked, using the Forward feature selection algorithm, and detection rates are plotted with increasing number of features (using the ranking provided by the Forward selection) on Fig. 6.

F5 case is different again: whatever the embedding rate used, steganalysis is surprisingly easy: using as low as 1 feature for all embedding rates, detection rate is 100% even with only 5% embedding rate.

The five other analyzed stego algorithms give rather different results:

• JPHS reaches a plateau in performance (within the standard deviation of 1%) for all embedding rates with 41 features and remains around that performance, even though performances are quite "unstable"

• OutGuess has this same plateau at 25 features and performances are not increasing anymore above that number of features (stll within the standard deviation of the results)

• StegHide can be considered to have reached the maximum result (within the standard deviation interval) at 60 features, even if for 5% embedding rate, performances keep on increasing while adding features. . .

• In the MM3 case, performances for embedding rates 10, 15 and 20% are very similar as are selected features. Performances stable at 40 features.

• Performances for MBSteg are stable using 70 features for embedding rates 15 and 20%. Only 30 are enough for embedding rate 5%. The case of embedding rate 10% has the classifier's performances increasing with the addition of features. Interestingly, the features retained by the forward selection for each embedding rate differ slightly, within one steganographic algorithm. Details about the features ranked as first by the Forward algorithm are discussed afterward.

C. Discussion

First, the global performances when using the reduced and sufficient feature sets mentioned in the results section above, are assessed. Note that feature selection for performing reverse-engineering of a steganographic algorithm is theoretically efficient only if the features are carrying different information (if two features represent the same information, the feature selection will select only one of them).

1) Reduced features sets: Based on the ranking of the features obtained by the Forward algorithm, it has been decided that once performances were within 1% of the best performance obtained (among all Forward tryouts for all different sets of features), the number of features obtained was retained as a "sufficient" feature set. Performances using reduced feature sets (proper to each algorithm and embedding rate) are first compared in Table I It should be noted that since the aim of the feature selection is to reduce as much as possible the feature set while keeping overall same performance, it is expected that within the standard deviation interval, the performance with the lowest possible number of features is behind the "maximum" one.

It remains possible, for the studied algorithms, as Fig. 6 shows, to find a higher number of features for which the performance is closer or equal to the maximum one -even though this is very disputable, considering the maximal 1% standard deviation interval when using 5000 images. But this is not the goal of the feature selection step of the methodology.

D. Feature sets analysis for reverse-engineering

Common feature sets have been selected according to the following rule: take the first common ten features (in the order ranked by the Forward algorithm) to each feature set obtained for each embedding rate (within one algorithm). It (f) StegHide Fig. 6. Performance in detection for all six stego algorithms versus the number of features, for the four embedding rates: black circles ( ) for 20%, green squares (2) for 15%, red crosses (×) for 10% and blue triangles ( ) for 5%. Features are ranked using the Forward selection algorithm. These plots are the result of a single run of the Forward algorithm. Plots do not have the same scale, vertically.

is hoped that through this selection, the obtained features will be generic regarding the embedding rate.

The F5 case is discarded, since only one feature makes it possible to achieve 99% of detection rate, that is the Global Histogram of coefficient 3: H(3) in the following notation.

Notations for the feature set used [START_REF] Pevny | Merging markov and dct features for multi-class jpeg steganalysis[END_REF] are first given for the original 23 features set, in Table II This set of 23 features is expanded up to a set of 193, by removing the L 1 norm used previously and keep all the values of the matrices and vectors. This results in the following 193 features set:

• A global histogram of 11 dimensions H(i), i = -5, 5

• 5 low frequency DCT histograms each of 11 dimensions h 21 (i) . . . h 31 (i), i = -5, 5

• 11 dual histograms each of 9 dimensions g -5 (i) . . . g 5 (i), i = 1, 9

• Variation of dimension 1 V • 2 blockinesses of dimension 1 B 1 , B 2 • Co-occurrence matrix of dimensions 5 × 5 C i,j , i = -2, 2 , j = -2, 2
Follows a discussion on the selected features for each steganographic algorithm.

Tables of selected feature sets are presented below, with an analysis of it for each algorithm. Fridrich's DCT features are not the only ones having a possible physical interpretation. They have been chosen here because it is believed that most of the features can be interpreted. The proposed short analysis of the weaknesses of stego algorithms is using this interpretation. 1) MM3: MM3 tends to be very sensitive to coefficients histograms features, which they do not preserve. Low DCT coefficients values (-1, +1) are found for MM3. Interestingly, co-occurrence coefficients react for MM3, as well as for JPHS (in the following) but only these two stego algorithms. 2) JPHS: JPHS seems not to preserve the low and medium frequencies coefficients and also not the frequency coherence (from the co-occurrence matrix). Also, DCT coefficients histograms are not preserved. 3) MBSteg: The features used (Table V) include global histograms with values 0,-2 and 2, which happens only because of the calibration in the feature extraction process. MBSteg preserves the coefficients' histograms but does not take into account a possible calibration. Hence, the unpreserved histograms are due to the calibration process in the feature extraction. Information leaks through the calibration process. and for low values (between -1 and +1).

h 21 (1) h 21 (2) h 12 (0) h 12 (2) h 22 (1) g -5 (1) C -2,-1 C -1,+1 C -1,+2 C +0,-2
h 12 (2) h 12 (3) h 13 (1) h 22 (1) h 31 (1) h 12 (2) C -2,-1 C -1,+1 C -1,+2 C +0,+1
From a general point of view, we can notice that all the analysed algorithms are very sensitive to statistics of low-pass calibrated DCT coefficients, which are represented by features h 21 and h 12 . This is not surprising since these coefficients contain a large part of the information of a natural image, their associated densities are likely to be modified by the embedding process. Note that fooling this process by embedding only inside the high-frequency coefficients is possible but reduces considerably the embedding payload.

V. CONCLUSIONS This paper has presented a methodology for the estimation of a sufficient number of images for a specific feature set using the standard deviation of the detection rate obtained by the classifier as a criterion (a Bootstrap technique is used for that purpose); the general methodology presented can nonetheless be extended and applied to other feature sets. The second step of the methodology aims at reducing the dimensionality of the data set, by selecting the most relevant features, according to a Forward selection algorithm; along with the positive effects of a lower dimensionality, analysis of the selected features is possible and gives insights on the steganographic algorithm studied. Three conclusions can be drawn from the methodology and experiments in this paper:

• Results on standard deviation for almost all studied steganographic algorithms have proved that the featurebased steganalysis is reliable and accurate only if a sufficient number of images is used for the actual training of the classifier used. Indeed, from most of the results obtained concerning standard deviation values (and therefore statistical stability of the results), it is rather irrelevant to possibly increase detection performance by 2% while working with a standard deviation for these same results of 2%.

• Through the second step of the methodology, the required number of features for steganalysis can be decreased. This with three main advantages: (a) performances remain the same if the reduced feature set is properly constructed; (b) the selected features from the reduced set are relevant and meaningful (the selected set can possibly vary, according to the feature selection technique used) and make reverse-engineering possible; (c) the weaknesses of the stego algorithm also appear from the selection; this can lead for example to improvements of the stego algorithm.

• The analysis on the reduced common feature sets obtained between embedding rates of the same stego algorithm, shows that the algorithms are sensitive to roughly the same features, as a basis. Meanwhile, when embedding rates get as low as 5%, or for very efficient algorithms, some very specific features appear.

Hence, the first step of the methodology is a requirement for any new stego algorithm, but also new feature sets/steganalyzers, willing to present its performances: a sufficient number of images for the stego algorithm and the steganalyzer used to test it has to be assessed in order to have stable results (that is, with a small enough standard deviation of its results to make the comparison with current state of the art techniques meaningful).

Also, from the second step of the methodology, the most relevant features can be obtained and make possible a further analysis of the stego algorithm considered, additionally to the detection rate obtained by the steganalyzer.

APPENDIX: FEATURES RANKED BY THE FORWARD ALGORITHM

As appendix are given the first 40 features obtained by the Forward ranking for each stego algorithm with 5% embedding rate. Only one embedding rate result is given for space reasons. 5% embedding rate results have been chosen since they tend to be different (in terms of ranked features by the Forward algorithm) from the other embedding rates and also because 5% embedding rate is a difficult challenge in terms of steganalysis; these features are meaningful for this kind of difficult steganalysis with these six algorithms. 5) H(4) g 2 (1) g 0 (1) g -3 (5) g 0 (9) g -3 (8) g -3 (3) g -5 (4) g -5 (5) C -2,-2 g -1 (6) g -2 (6)

g 4 (3) C -1,-1 C -1,0 g -2 (7) C -1,1

Fig. 2 .

 2 Fig.2. Scheme of the possible reverse-engineering on an unknown stego algorithm, by using feature selection for identification of the specific weaknesses.

Fig. 4 .

 4 Fig. 4. Structure of a classical Single Layer Feedforward Neural Network (SLFN). The input values (the data) X = (x 1 , . . . , x N ) are weighted by the W coefficients. A possible bias B (not on the figure) can be added to the weighted inputs w i x i . An activation function f taking this weighted inputs (plus bias) as input is finally weighted by output coefficients β to obtain the output Y = (y 1 , . . . , y N ).

  .

		5%	#	10%	#
	F5	100.0	1	100.0	1
	JPHS	90.7	41	92.1	21
	MBSteg	63.3	57	70.9	93
	MM3	78.00	81	86.2	49
	OutGuess	81.2	65	93.2	49
	Steghide	82.3	149	91.2	89
		15%	#	20%	#
	F5	100.0	1	100.0	1
	JPHS	93.7	41	97.3	25
	MBSteg	83.5	73	88.5	69
	MM3	86.6	57	86.6	73
	OutGuess	98.8	33	100.0	29
	Steghide	96.4	73	99	73

TABLE I PERFORMANCES

 I FOR OP-ELM LOO FOR THE BEST FEATURES SET ALONG WITH THE SIZE OF THE REDUCED FEATURE SET (#). PERFORMANCES USING THE REDUCED SET ARE WITHIN THE 1% RANGE OF STANDARD DEVIATION OF THE BEST RESULTS. THE SIZE OF THE SET HAS BE DETERMINED TO BE THE SMALLEST POSSIBLE ONE GIVING THISPERFORMANCE.

  : 21 /||h 21 ||,h 12 /||h 12 ||,h 13 /||h 13 ||, for 5 DCT Modes h 22 /||h 22 ||,h 31 /||h 31 || Dual histogram for g -5 /||g -5 ||,g -4 /||g -4 ||,. . . , 11 DCT values g 4 /||g 4 ||,g 5 /||g 5 ||

	Functional/Feature	Functional F
	Global histogram	H/||H||
	Individual histogram h Variation V
	L1 and L2 blockiness	B 1 , B 2
	Co-occurrence	N 00 , N 01 , N 11

TABLE IV COMMON

 IV FEATURE SET FOR JPHS.

  Extreme values histograms are mostly used (values -2, -1) in the feature set for OutGuess (TableVI) and a clear weak point. The calibration process has indeed been of importance since the histogram of value 0 has been taken into account.

		TABLE VI
	COMMON FEATURE SET FOR OUTGUESS.
	4) OutGuess: h 21 (-1) h 21 (0) h 12 (1)	h 12 (3)	h 13 (-1)
	h 13 (0)	h 13 (1) h 22 (1)	h 31 (-1)	h 31 (0)

  TABLE VIII THE 40 FIRST FEATURES RANKED BY THE FORWARD ALGORITHM FOR THE F5 ALGORITHM AT 5% EMBEDDING RATE.

	H(3)	g 5 (1)	H(-3)	g -2 (7)	h 13 (0)	C -2,-1	H(4)	h 12 (-3)	g 2 (4)	g -1 (4)
	g 3 (1)	g0(3)	h 21 (-3)	g -4 (2)	g 5 (6)	g -3 (2)	g 2 (7)	h 13 (-3)	H(-1)	H(-4)
	H(-2)	g 5 (2)	H(1)	h 12 (-2)	h 31 (-1)	h 12 (5)	h 12 (-4)	h 13 (-2)	H(5)	H(2)
	h 21 (5)	h 12 (4)	h 13 (-5)	g -2 (1)	h 21 (0)	h 21 (-5)	g 1 (3)	h 21 (3)	h 22 (-1)	h 21 (2)
	g 0 (4)	h 22 (0)	C 1,0	B1	H(1)	h 21 (0)	g 1 (4)	g 0 (8)	g -2 (9)	g -2 (5)
	g 4 (5)	g 0 (5)	g 1 (9)	g -1 (2)	B2	g 2 (8)	C 0,0	h 31 (5)	g 0 (9)	h 22 (1)
	g -2 (2)	g -1 (7)	g -3 (8)	g 0 (1)	h 31 (-3) h 21 (-1)	h 22 (-1)	g -4 (6) C -1,-2	g 5 (7)
	h 12 (-5)	g -5 (8)	h 21 (2)	g 0 (7)	h 12 (-2) h 22 (-4)	h 31 (0)	C 0,2	H(2)	g 5 (5)

TABLE IX THE

 IX 40 FIRST FEATURES RANKED BY THE FORWARD ALGORITHM FOR THE JPHS ALGORITHM AT 5% EMBEDDING RATE.

	g -2 (1)	H(2)	g -4 (7)	h 13 (1)	h 22 (1)	C 2,-2	C -1,-1	h 31 (1)	g 4 (7)	g -2 (4)
	h 21 (0)	h 31 (-4) h 21 (-4)	C 0,2	C 1 2	h 31 (-1)	H(0)	h 21 (3)	g -5 (6)	h 22 (-3)
	h 13 (-1)	C 2,0	C 1,2	g 5 (6)	C -2,-1	g -3 (6)	g 5 (4)	g -2 (7)	g -1 (7)	g -4 (8)
	h 22 (-1)	g 2 (1)	g 0 (8)	h 22 (-5)	H(-2)	h 12 (-4)	g 5 (5)	h 12 (-2)	g 2 (4)	h 21 (-3)

TABLE X THE

 X 40 FIRST FEATURES RANKED BY THE FORWARD ALGORITHM FOR THE MBSTEG ALGORITHM AT 5% EMBEDDING RATE.

	C -1,-1	h 13 (-1)	C 0,-2	C 1,1	g 0 (9)	C 2,0	h 21 (-1)	h 13 (1)	g -3 (2)	C 1 0
	H(-2)	g 4 (4)	g 2 (2)	C -2,0	C 0,-1	C -1,-2	g -2 (3)	h 22 (-3)	g 2 (3)	h 13 (3)
	h 31 (-1)	g -1 (9)	g -2 (8)	g 0 (7)	h 21 (-5)	h 21 (3)	C -1,1	g -1 (3)	g 5 (3)	h 31 (1)
	g 0 (3)	B1	C -2,1	B2	g -4 (6)	C 0,2	H(-1)	g 2 (5)	h 13 (0)	g 2 (7)

TABLE XI THE

 XI 40 FIRST FEATURES RANKED BY THE FORWARD ALGORITHM FOR THE MM3 ALGORITHM AT 5% EMBEDDING RATE.

	h 13 (0)	C 0,-1	C -2,0	H(-2)	B1	C 0,-2	g 0 (7)	h 31 (-3)	C -2,-1	g 0 (2)
	B2	H(-1)	g -2 (2)	h 13 (-1) h 22 (-1)	h 22 (0)	h 12 (-3)	g -2 (5)	g 1 (8)	h 21 (-2)
	g -2 (9)	g 1 (1)	H(

TABLE XII THE

 XII 40 FIRST FEATURES RANKED BY THE FORWARD ALGORITHM FOR THE OUTGUESS ALGORITHM AT 5% EMBEDDING RATE.TABLE XIII THE 40 FIRST FEATURES RANKED BY THE FORWARD ALGORITHM FOR THE STEGHIDE ALGORITHM AT 5% EMBEDDING RATE.

	C 0,-1	g 0 (2)	C 0,2	C 2,-2	B1	B2	C 1,1	C 0,-2	C -2,2	h 13 (-1)
	g -5 (3)	h 21 (-3)	C 0,1	h 13 (0)	C 1,-1	h 31 (-1) g -3 (3)	g 3 (6)	h 31 (-2)	g 1 (3)
	h 22 (1)	C -2,-2	g -4 (4)	h 13 (1)	C -2,0	g 1 (4)	C 2,1	H(-1)	C 2,2	h 22 (5)
	g 2 (5)	C -1,-1	g 1 (9)	C 2,0	g 2 (7)	g -1 (1)	h 31 (5)	H(-2)	h 21 (1)	g -2 (9)

ACKNOWLEDGEMENTS

The authors would like to thank Jan Kodovsky and Jessica Fridrich for their implementation of the DCT Feature Extraction software.


Xavier Guyon
email: guyon@univ-paris1.fr

Cécile Hardouin
email: hardouin@univ-paris1.fr

Forward recursions and normalizing constant for Gibbs fields

Keywords: Gibbs distribution, interaction potential, Markov Chain, Markov field, marginal law, normalizing constant, Ising model

Maximum likelihood parameter estimation is frequently replaced by various techniques because of its intractable normalizing constant. In the same way, the literature displays various alternatives for distributions involving such unreachable constants.

In this paper, we consider a Gibbs distribution π and present a recurrence formula allowing a recursive calculus of the marginals of π and in the same time its normalizing constant.

The numerical performance of this algorithm is evaluated for several examples, particularly for an Ising model on a lattice.

Introduction

Usually, the normalizing constant of a discrete probability π distribution involves high dimensional summation such that the direct evaluation of these sums becomes quickly infeasible in practice. For example, the direct calculation of the normalizing constant of the Ising model on a 10×10 grid involves summation over 2 100 terms. The problem is bypassed for instance by evicting the distribution of interest by an alternative as, for example in spatial statistics, replacing the likelihood for the conditional pseudo likelihood (see [START_REF] Besag | Spatial interactions and the statistical analysis of lattice systems[END_REF]). Another solution consists of estimating the normalizing constant; a number of techniques have been proposed for this approximate evaluation, see for example [START_REF] Pettitt | Efficient calculation of the normalizing constant of the autologistic and related models on the cylinder and lattice[END_REF] and [START_REF] Moeller | An efficient Markov chain Monte Carlo method for distributions with intractable normalizing constants[END_REF] for efficient Monte Carlo methods. [START_REF] Bartolucci | A recursive algorithm for Markov random fields[END_REF] present a recursive algorithm computing the likelihood of a Markov random field in the form of a product of conditional probabilities. Reeves (2004) propose efficient computation of the normalizing constant for a factorisable model, that is when π can be written as a product of factors.

In this paper we give more specific results when we specify π as a Gibbs distribution. We derive some results of Khaled (see [START_REF] Khaled | A multivariate generalization of a Markov switching model[END_REF], Khaled 2 (2008)) giving an original linear recursive calculation of the marginal distributions in the case of a particular distribution π of Z = (Z 1 , Z 2 , • • • , Z T ) used in econometrics for the modelling of a hidden regime (see [START_REF] Lovinson | A matrix-valued Bernoulli distribution[END_REF]). An interesting consequence is to ease the calculation of the normalizing constant of π. We generalize Khaled's results, noticing that if π is a Gibbs distribution on T = {1, 2, • • • , T }, therefore π is a Markov field on T and it is easy to manipulate its conditional distributions. This approach allows to extend the recurrences given by Khaled to general Gibbs distribution π (either spatial or temporal) on a general finite state space E; those recursions yield exact calculation of the marginal distributions

π t of (Z 1 , Z 2 , • • • , Z t ), 1 ≤ t ≤ T as well as the normalization constant C of π.
First, we recall in section 1 some basic properties about Gibbs and Markov fields for

Z = (Z 1 , Z 2 , • • • , Z T ), a sequence of joint distribution π on E T where E is a finite state space.
The main result of this paper is presented in section 2, where we give forward recursion for the marginal distributions and the application to the calculation of the normalization constant. We present some simple examples and give the computing times for carrying out the normalizing constant. In section 3, we extend the results to general Gibbs fields, in the sense of temporal Gibbs fields as well as spatial Gibbs fields. Finally we evaluate the numerical performance of our method to compute the normalizing constant for a spatial Ising model on a lattice m × T .

Markov chain and Markov field properties for a Gibbs field

Let T > 0 be a fix positive integer,

E = {e 1 , e 2 , • • • , e N } a finite state space with N elements, Z = (Z 1 , Z 2 , • • • , Z T ) a temporal sequence with joint distribution π on E T .
We assume that π is a Gibbs distribution with an energy U T described by singletons potentials (θ s ) s=1,T and pairs potentials (Ψ s ) s=2,T , that is, denoting

z(t) = (z 1 , z 2 , • • • , z t ) for 1 ≤ t ≤ T , π(z(T )) = C exp U T (z(T )) with C -1 = z(T )∈E T exp U T (z(T )) where (2•1) U t (z(t)) = s=1,t θ s (z s ) + s=2,t Ψ s (z s-1 , z s ) for 2 ≤ t ≤ T , and U 1 (z 1 ) = θ 1 (z 1 ).

2•1 Markov field property

The neighbourhood graph on T = {1, 2, • • • , T } associated to π is the so-called 2 nearest neighbours system. Moreover, we consider π as a bilateral Markov random field equipped with the 2 nearest neighbours system (see [START_REF] Kindermann | Markov random fields and their applications[END_REF], [START_REF] Guyon | Random Fields on a Network: Modeling, Statistics, and Applications[END_REF]

): if 1 < t < T , π(z t | z s , 1 ≤ s ≤ T and s = t) = π(z t | z t-1 , z t+1 ) (2•2) Indeed: π(z t | z s , 1 ≤ s ≤ T and s = t) = π(z 1 , z 2 , • • • , z T ) u∈E π(z 1 , z 2 , • • • , z t-1 , u, z t+1 , ..z T ) = exp{θ t (z t ) + Ψ t (z t-1 , z t ) + Ψ t+1 (z t , z t+1 )} u∈E exp{θ t (u) + Ψ t (z t-1 , u) + Ψ t+1 (u, z t+1 )} = π(z t | z t-1 , z t+1 ).
Therefore, the non causal conditional distribution π(z t | z t-1 , z t+1 ) can be easily computed as soon as N , the cardinal of E, remains rather small.

2•2 Markov chain property

Proposition 1 Z is a Markov chain: π(z t | z s , s ≤ t -1) = π(z t | z t-1 ) if 1 < t ≤ T . Proof : π(z t | z s , s ≤ t -1) = πt(z 1 ,z 2 ,••• ,zt) π t-1 (z 1 ,z 2 ,••• ,z t-1 ) . Let's identify π t (z 1 , z 2 , • • • , z t )
. For 1 ≤ s < t ≤ T , and using the notation u t s = (u s , u s+1 , • • • , u t ), we write:

π t (z 1 , z 2 , • • • , z t ) = u T t+1 ∈E T -t π(z 1 , z 2 , • • • , z t , u T t+1 ) (2•3) = C exp{ s=1,t θ s (z s ) + s=2,t Ψ s (z s-1 , z s )} × exp{θ * t (z t )} where exp{θ * t (z t )} = u T t+1 exp{ T s=t+1 θ s (u s ) + Ψ t+1 (z t , u t+1 ) + T t+2 Ψ s (u s-1 , u s )}. (2•4)
We deduce that the conditional (to the past) distribution of z t is given by:

π(z t | z s , s ≤ t -1) = exp{θ t (z t ) + θ * t (z t ) -θ * t-1 (z t-1 ) + Ψ t (z t-1 , z t )} = π(z t | z t-1 ). (2•5)
Remarks:

1. Following the proof above, we write again the marginal distribution

π t (z 1 , z 2 , • • • , z t ) as π t (z(t)) = C exp{ s=1,t-1 θ s (z s ) + [θ t (z t ) + θ * t (z t )] + s=2,t Ψ s (z s-1 , z s )}.
Therefore, the marginal field (Z 1 , • • • , Z t ) is also a Gibbs field with the 2 nearest neighbours system, associated to the same potentials, except for the last singleton potential which equals θ t (z t ) = θ t (z t ) + θ * t (z t ).

2. An important difference appears between formula (2•5) and (2•2): indeed, (2•2) is computationally feasible, when (2•5) is not, because of the summation over u T t+1 which is of complexity N T -t .

3 Recursion over marginal distributions

3•1 Future-conditional contribution Γ t (z(t)) For t ≤ T -1, the distribution conditionally to the future π(z 1 , z 2 , • • • , z t | z t+1 , z t+2 , • • • , z T ), depends only on z t+1 : π(z 1 , z 2 , • • • , z t | z t+1 , z t+2 , • • • , z T ) = π(z 1 , z 2 , • • • , z T ) u t 1 ∈E t π(u t 1 , z t+1 , ..z T ) = exp {U t (z 1 , ..., z t ) + Ψ t+1 (z t , z t+1 )} u t 1 ∈E t exp s=1,t θ s (u s ) + Ψ t+1 (u t , z t+1 ) + s=2,t Ψ s (u s-1 , u s ) = π(z 1 , z 2 , • • • , z t | z t+1 ).
We can write this conditional distribution on the following feature:

π(z 1 , z 2 , • • • , z t | z t+1 ) = exp {U t (z 1 , ..., z t ) + Ψ t+1 (z t , z t+1 )} u t 1 ∈E t exp {U t (u 1 , ..., u t ) + Ψ t+1 (u t , z t+1 )} , i.e. π(z 1 , z 2 , • • • , z t | z t+1 ) = C t (z t+1 ) exp U * t (z 1 , z 2 , • • • , z t ; z t+1 )
where U * t is the future-conditional energy defined by:

U * t (z 1 , z 2 , • • • , z t ; z t+1 ) = U t (z 1 , z 2 , • • • , z t ) + Ψ t+1 (z t , z t+1 ), (3•1)
and

C t+1 (z t+1 ) -1 = u t 1 ∈E t exp {U * t (u 1 , ..., u t ; z t+1 )}. Then, for i = 1, N : π(z 1 , z 2 , • • • , z t | z t+1 = e i ) = C t (e i )γ t (z 1 , z 2 , • • • , z t ; e i ) where γ t (z(t); e i ) = exp U * t (z(t); e i ).
Thus, we define γ t (z(t); e i ) as the contribution to the distribution π t of Z(t) conditionally to the future z t+1 = e i .

Definition 1 For t ≤ T -1, the vector Γ t (z(t)) of the future-conditional contributions is defined by the vector of R N with i-th component, 1 ≤ i ≤ N :

(Γ t (z(t))) i = γ t (z(t); e i ).
For t = T , there is no conditional future and Γ T (z(T )) is the constant vector of components γ T (z(T )) = exp U T (z(T )). The definition of Γ T (z(T )) is analogous to the one of Γ t for t ≤ T -1 with the convention Ψ T +1 ≡ 0. Still with this convention, we define for 1 ≤ t ≤ T the matrix A t of size N × N with general term:

A t (i, j) = exp{θ t (e j ) + Ψ t+1 (e j , e i )}, et i, j = 1, N. (3•2)
Let us note that A T has constant columns A T (i, j) = exp θ T (e j ) for i, j = 1 to N . Then we get the fundamental recurrence:

Proposition 2 For all 2 ≤ t ≤ T , z(t) = (z 1 , z 2 , • • • , z t ) ∈ E
t and e i ∈ E, we have:

γ t (z(t -1), e j ; e i ) = A t (i, j) × γ t-1 (z(t -1); e j ) , (3•3) and z t ∈E Γ t (z(t -1), z t ) = A t Γ t-1 (z(t -1)). (3•4) Proof : (i) Let us consider the case 2 ≤ t ≤ T -1. The energy U t verifies U t (z(t -1), z t ) = U t-1 (z(t -1)) + θ t (z t ) + Ψ t (z t-1 , z t ) ; therefore (3•1) ensures for all (z t , z t+1 ) = (a, b) ∈ E 2 : U * t (z(t -1), a; b) = U t-1 (z(t -1)) + θ t (a) + Ψ t (z t-1 , a) + Ψ t+1 (a, b) = U * t-1 (z(t -1); a) + {θ t (a) + Ψ t+1 (a, b)}.
This implies the recurrence (3•3), with (e j , e i ) = (a, b) = (z t , z t+1 ). The summation over z t = e j gives the component i on the left hand side of (3•4). That ensures the result.

(ii

) For t = T , since Ψ T +1 ≡ 0, then U * T (z(T -1), a) = U * T -1 (z(T -1); a) + θ T (a). Therefore we obtain (3•4) for A T (i, j) = exp θ T (e j ).

3•2 Forward recursion for marginal distributions and normalization constant

Let us define the following row vectors 1 × N : E 1 = B T = (1, 0, • • • , 0), and (B t ) t=T,2 the sequence specified by the forward recursion:

B t-1 = B t A t if t ≤ T.
We also denote

K 1 = z 1 ∈E Γ 1 (z 1 ) ∈ R N . The evaluation of K 1 is easy if N is not too large, since its i-th component is K 1i = z 1 ∈E exp{θ 1 (z 1 ) + Ψ 2 (z 1 , e i )}.
We give below the main result of this paper.

Proposition 3 Marginal distributions π t and calculation of the normalization constant C.

(1) For 1 ≤ t ≤ T :

π t (z(t)) = C × B t Γ t (z(t)). (3•5)
(2) The normalization constant C of the joint distribution π verifies:

C -1 = E 1 A T A T -1 • • • A 2 K 1 . (3•6)
Proof :

(1) Let us prove (3•5) by descending recurrence. For t = T , the equality is verified since,

π(z 1 , z 2 , • • • , z T ) = π T (z(T )) = C exp U T (z(T )) = C × E 1 Γ T (z(T )).
Let us assume that (3•5) is verified for t, 2 ≤ t ≤ T . We use (3•4) which gives:

π t-1 (z(t -1)) = z t ∈E π t (z(t -1), z t ) = C × B t { z t ∈E Γ t (z(t -1), z t )} = C × B t A t Γ t-1 (z(t -1)) = C × B t-1 Γ t-1 (z(t -1)). (2) Refering to (3•5), we have π 1 (z 1 ) = C ×B 1 Γ 1 (z 1 ). The summation over z 1 gives C × B 1 K 1 = 1. The result ensures from the equality B 1 = E 1 A T A T -1 • • • A 2 .

Remarks :

1 -The formula (3•6) allowing the calculus of C reduces in the case of time invariant potentials. In that case, we have A t ≡ A for 1 ≤ t ≤ T -1, with A(i, j) = exp{θ(e j ) + Ψ(e j , e i )}, A T (i, j) = exp θ(e j ) and

C -1 = E 1 A T A T -2 K 1 . (3•7)
Therefore, if the size N of E allows the diagonalization of the matrix A, we can achieve the calculus of the constant C independently of the temporal dimension T .

We give below two examples to illustrate our method. In each example, we consider Gibbs models with time independent potentials. For increasing values of T and various methods, we compare the times necessary for the computation of the normalizing constant C -1 . When possible, we also give the computing times for the ordinary summation method.

In all examples, we have used the Matlab software and times are given in seconds.

Example 1 : binary temporal model The state space is E = {0, 1} = {e 1 , e 2 } with N = 2 states, θ t (z t ) = αz t and Ψ t+1 (z t , z t+1 ) = βz t z t+1 for t ≤ T -1. We have :

A t = A = 1 e α 1 e α+β for t = 1, T -1, A T = 1 e α 1 e α , E 1 = 1 0 , K 1 = 1 + e α 1 + e α+β .
We present in Table 1 the times for the computation of C -1 for various values of T and the following methods : (1) calculation of

C -1 = E 1 A T A T -2 K 1 ; (2) C -1
is obtained by direct summation on E T using a simple loop, or using a bitmap dodge which computes simultaneously the 2 T elements of E T (2 bis). We stopped computing C -1 by summation for T > 25.

Meth 1 : Computing times of C -1 for the binary temporal Gibbs distribution with parameters α = 1, β = -0.8.

We see that the computing times of C -1 = E 1 A T A T -2 K 1 are negligible for T < 700 while the summation method becomes quickly infeasible.

Example 2 : bivariate binary temporal model E = {0, 1} 2 (N = 4 states), and Z(T ) is the anisotropic Ising model:

π 2 ((x 1 , y 1 ), (x 2 , y 2 )) = C exp{αx 1 + βy 1 + γx 1 y 1 + αx 2 + βy 2 + γx 2 y 2 + δ(x 1 x 2 + y 1 y 2 )}. E 1 A T = A T (1, •) is A T 's
first row with A T (1, j) = exp{αx j + βy j + γx j y j }; A is of size 4 × 4 defined by A(i, j) = exp{αx j + βy j + γx j y j + δ(x i x j + y i y j )}, and K 1 's i-th component equals z=(x,y)∈E exp{αx + βy + γxy + δ(xx i + yy i )}. We computed C -1 in two ways, first using the the powers of A i.e.

C -1 = E 1 A T A T -2 K 1 , then using its diagonalization C -1 = E 1 A T P D T -2 P -1 K 1 .
We took parameters α = 1, β = -0.8, γ = -0.5, δ = 0.04. We were able to calculate C -1 = 9.9491e + 006 for T = 430 and then stopped for larger T since the software treats C -1 as equals infinity. The computing times are null for both methods, which means that we get instantaneously C -1 's value, and the size of A is still too small to distinguish computations using power or diagonalization of A.

General Gibbs fields

4•1 Temporal Gibbs fields

The previous results can be extended to general temporal Gibbs models. As an illustration, let us consider the following model, characterized by the energy:

U T (z(T )) = s=1,T θ s (z s ) + s=2,T Ψ 1,s (z s-1 , z s ) + s=3,T Ψ 2,s (z s-2 , z s ).
The joint distribution π defines a bilateral Markov field with the 4 nearest neighbours system, and conditional distributions:

π(z t | z s , 1 ≤ s ≤ T et s = t) = = exp{θ t (z t ) + Ψ 1,t (z t-1 , z t ) + Ψ 1,t+1 (z t , z t+1 ) + Ψ 2,t (z t-2 , z t ) + Ψ 2,t+2 (z t , z t+2 )} u∈E exp{θ t (u) + Ψ 1,t (z t-1 , u) + Ψ 1,t+1 (u, z t+1 ) + Ψ 2,t (z t-2 , u) + Ψ 2,t+2 (u, z t+2 )} = π(z t | z t-1 , z t+1 , z t-2 , z t+2 ).
Z is also a Markov chain of order 2 with:

π(z t | z s , s ≤ t -1) = exp{θ t (z t ) + θ * t (z t ) -θ * t-1 (z t-1 ) + Ψ 1,t (z t-1 , z t ) + Ψ 2,t (z t-2 , z t ) + θ * * t (z t-1 , z t ) -θ * * t-1 (z t-2 , z t-1 )} = π(z t | z t-1 , z t-2 )
where θ * t (z t ) is given by (2•4) and

θ * * t (z t-1 ) = u T t+1 exp{Ψ 2,t+1 (z t-1, u t+1 ) + Ψ 2,t+2 (z t , u t+2 ) + T t+3 Ψ 2,s (u s-2, u s )}. For t ≤ T -2, the conditional distribution π(z(t) | z t+1 , z t+2 , • • • , z T ) depends only on (z t+1 , z t+2 ): π(z(t) | z t+1 , z t+2 ) = C t (z t+1 , z t+2 ) exp U * t (z(t); z t+1 , z t+2 ),
where we name U * t (z(t); z t+1 , z t+2 ) the future conditional energy:

U * t (z(t); z t+1 , z t+2 ) = U t (z(t)) + Ψ 1,t+1 (z t , z t+1 ) + Ψ 2,t+1 (z t-1 , z t+1 ) + Ψ 2,t+2 (z t , z t+2 ), and C t (z t+1 , z t+2 ) -1 = u t 1 ∈E t exp U * t (u 1 , .., u t ; z t+1 , z t+2 ).
For a, b and c ∈ E, it is easy to verify that:

U * t (z(t -1), a; (b, c)) = U * t-1 (z(t -1); (a, b)) + θ t (a) + Ψ 1,t+1 (a, b) + Ψ 2,t+2 (a, c).
With the convention Ψ 1,s ≡ Ψ 2,s ≡ 0 for s > T , we define:

• for t ≤ T -2, the vector Γ t (z(t)) of the conditional contributions, conditionally to the future (z t+1 , z t+2 ) = (e i , e j ), i, j = 1, N, by the N 2 × 1 vector of components (i, j): γ t (z(t); (e i , e j )) = exp U * t (z(t); e i , e j );

• Γ T -1 (z(T -1)
), the vector of the contributions conditionally to the future z T = e i :

(Γ T -1 (z(T -1)) i = exp{U T -1 (z(T -1)) + Ψ 1,T (z T -1 , e i ) + Ψ 2,T (z T -1 , e i )};
• Γ T (z(T )) the constant vector of components exp{U T (z(T )}.

In the same way, we define the matrix A of size N 2 × N 2 whose non zero components are:

A t ((i, j), (k, i)) = exp{θ t (e k ) + Ψ 1,t+1 (e k , e i ) + Ψ 2,t+2 (e k , e j )}
Like in section 3•2, we obtain a recurrence formula on the γ t :

γ t (z(t -1), e k ; (e i , e j )) = A t ((i, j), (k, i)) × γ t-1 (z(t -1); (e k , e i ))
together with the statement (3•4) on the contributions Γ t (z(t)). In this context, (Z t , t = 1, T ) is a Markovian process with memory 2 and Y t = (Z t , Z t+1 ) a bivariate Markov chain for which we get the results (3•5) and (3•6).

4•2 Spatial Gibbs fields

Let us consider Z t = (Z (t,i) , i ∈ I), where I = {1, 2, • • • , m}, and

Z (t,i) ∈ F (Z t ∈ E = F m ). Then Z = (Z s , s = (t, i) ∈ S
) is a spatial field on S = T × I. We note again z t = (z (t,i) , i ∈ I), z(t) = (z 1 , .., z t ) and z = z(T ).

Without loss of generality, we suppose that the distribution π of Z is a Gibbs distribution with translation invariant potentials

Φ A k (•), k = 1, K associated to a family of subsets {A k , k = 1, K} of S, Φ A k (z) depending only on z A k , the layout of z over A k .
Then π is characterized by the energy:

U (z) = k=1,K s∈S(k) Φ A k +s (z), with S(k) = {s ∈ S s.t. A k + s ⊆ S}.
For A ⊆ S, we define the height of A by H(A) = sup{|u -v| , ∃(u, i) and (v, j) ∈ A}, and H = sup{H(A k ), k = 1, K} the biggest height of the potentials. With this notation, we write the energy U as the following:

U (z) = H h=0 T t=h+1 Ψ(z t-h , • • • , z t ) with Ψ(z t-h , • • • , z t ) = k:H(A k )=h s∈S t (k) Φ A k +s (z) where S t (k) = {s = (u, i) : A k + s ⊆ S and t -H(A k ) ≤ u ≤ t}.
(Z t ) is a Markov field with the 2H-nearest neighbours system but also a Markov process with memory H; Y t = (Z t-H , Z t-H+1 , • • • Z t ), t > H, is a Markov chain on E H for which we get the results (3•5) and (3•6).

4•3 Computing the normalization constant for the Ising model

We specify here the calculus of C -1 in the case of a translation invariant potentials Ising model [START_REF] Kindermann | Markov random fields and their applications[END_REF], [START_REF] Guyon | Random Fields on a Network: Modeling, Statistics, and Applications[END_REF]

). Let S = T × I ={1, 2, • • • , T } × {1, 2, • • • ,
m} be the set of sites, and F = {-1, +1} the state space. We consider Z = (Z (t,i) , (t, i) ∈ S) a Markov field on S with the four nearest neighbours system, a site (t, i) being a neighbour of (s, j) if (t, i) -(s, j) 1 = 1. The joint distribution π of Z is characterized by the singletons and pairs potentials:

Φ t,i (z) = α z (t,i) for (t, i) ∈ S, Φ {(t,i),(t,i+1)} (z) = β z (t,i) z (t,i+1) for 1 ≤ i ≤ m -1, and Φ {(t-1,i),(t,i)} (z) = δ z (t-1,i) z (t,i) for 2 ≤ t ≤ T.
Forgetting the spatial dimension, we consider the state z t = (z (t,i) , i = 1, m) ∈ E = {-1, +1} m . Then Z is a temporal Gibbs field with the following translation invariant potentials:

θ t (z t ) = θ(z t ) = α i=1,m z (t,i) + β i=1,m-1 z (t,i) z (t,i+1) , Ψ t (z t-1 , z t ) = Ψ(z t-1 , z t ) = δ i=1,m z (t-1,i) z (t,i) , 2 ≤ t ≤ T.
Let us give some notations associated to c = (c i , i = 1, m) and d = (

d i , i = 1, m), two states of E = {-1, +1} m ; we first introduce n + (c) = {i ∈ I : c i = +1} and n -(c) = {i ∈ I : c i = -1} (n + (c) + n -(c) = m), then v + (c) = {i = 1, m -1 : c i = c i+1 } and v -(c) = {i = 1, m -1 : c i = c i+1 } (v + (c) + v + (c) = m -1), and finally n + (c, d) = {i ∈ I : c i = d i }, n -(c, d) = {i ∈ I : c i = d i } (n + (c, d) + n + (c, d) = m).
Since the potentials are invariant, the matrix A t given by (3•2) does not depend on t, we have

A t = A, t ≤ T -1, the 2 m × 2 m matrix whose general term is, A(a, b) = exp{α(n + (b) -n -(b)) + β(v + (b) -v -(b)) + δ(n + (a, b) -n -(a, b))}, a, b ∈ E. Moreover, for t = T, A T (a, b) = exp{α(n + (b) -n -(b)) + β(v + (b) -v -(b))}. Therefore, the normalization constant of π is C -1 = E 1 A T A T -2 K 1 .
Since K 1 is given by a summation over E = {-1, +1} m , and the size A is 2 m × 2 m , this formula is practically useful for m not too big.

Example 3 Table 2 gives computing times for the normalizing constant for the Ising model above with parameters α = 0.15, β = 0.05, δ = -0.08. In the case T = m = 2, the model is the one given in example 2 with F = {0, 1}. We consider varying values of T and m = 10, that is we work with vectors and square matrices of size 2 10 but without theoretical constraints on the size of T = {1, 2, • • • , T }. We compute C -1 using the powers of matrix A or its diagonalization. We observe that it's computationally more efficient to compute the powers A T -2 rather than to use the diagonalization of A. Indeed, the diagonalization procedure itself is expensive for large size matrices.

m = 10 C -1 = E 1 A T A T -2 K 1 C -1 = E 1 A T P D T -2 P -1 K 1 Value C -1 T = 2 0.

4•4 Other generalizations

First, we presented the forward recursion for Gibbs distributions with singletons and pairs potentials. The results can be extended to larger potentials (triple or more).

Another extension is to consider variable state spaces. The recurrence (3•4) and properties (3•5), (3•6) hold for different state spaces E t of the components Z t ; in this case, the associated matrices A t involved in (3•6)are not necessarly square.

Finally, we can extend the results to embedded sets T ; for instance let us consider the decreasing sequence

T = S Q ⊃ S Q-1 ⊃ • • • ⊃ S 1 of parts of T = {1, 2, • • • , T };
similarly to the former future-conditional contributions, we define the contributions γ q (z(S q ); z(T S q )), conditionally to the outer layout z(T S q ). Let us give the following example: we assume S q = S q-1 ∪ ∂S q-1 for q = 1, Q -1. Then, we obtain the conditional energy U * q (z(S q ); z(∂S q )) = U * q-1 (z(S q-1 ); z(∂S q-1 )) + ∆ q (z(∂S q-1 )); z(∂S q )) with ∆ q (z(∂S q-1 )); z(∂S q )) = u∈∂S q-1 θ u (z u )+ u∈∂S q-1 ,v∈∂Sq,<u,v> Ψ {u,v} (z u , z v ). then we define the matrices A q by A q (∂S q ; ∂S q-1 ) = exp ∆ q (z(∂S q-1 ); z(∂S q )).

As an illustration let us set the following decreasing sequence

T = S T -1 = {1, 2, • • • , T }, S T -2 = {1, 2, • • • , T -1}, • • • , S 2 = {1, 2, 3} and S 1 = {2}. For q = T -1, • • • , 3,
the conditional contributions and the matrices A are defined in the usual way (3•1) and (3•2), while for q = 2, A 2 is a N × N 2 matrix with A 2 (z 4 , (z 1 , z 3 )) = exp{θ 1 (z 1 ) + θ 3 (z 3 ) + Ψ 4 (z 3 , z 4 )}.

Conclusion

The technique proposed in this paper to evaluate marginals and normalisation constant is applicable to Markov chains, Markov fields. It overcomes the need to resort to approximate alternatives when one wants to evaluate the likelihood or the normalizing constant of a Gibbs field. It makes feasible the exact evaluation of the normalizing constant for moderate set of sites, eliminating any Monte Carlo procedure, variational scheme, or (approximate) Bayesian computations steps. As a statistical consequence, it is therefore possible to perform true maximum likelihood estimation for such Gibbs field. Another application is to allow exact simulation.

We gave several illustrations of computing times for the normalizing constant. For one dimensional two states Gibbs fields, we are able to compute instantaneously the normalizing constant for a sequence of length 700, as well as for a sequence with four states and of length 400. We could keep computing for bigger lengths using another software. However, our main goal here is to provide a new method and we let the users choose their way of programming.

For spatial processes, we have computed the normalizing constant for an Ising model on a lattice 10 × 100 in 20 seconds. Following the discussion above, we could increase T, one of the side of the lattice; the limitation of the procedure ensures from the manipulation of N m × N m matrices. So the method seems to fail for large square lattices. As a comparison, [START_REF] Pettitt | Efficient calculation of the normalizing constant of the autologistic and related models on the cylinder and lattice[END_REF] compute the normalizing constant for an autologistic model defined on a cylinder lattice for which the smallest row or column is not greater than 10. They suggest to split a large lattice into smaller sublattices along the smallest row or column. A similar idea could apply here.

Table 2 :

 2 Computing times of C -1 for an Ising model on a lattice 10 × T for various T.

		3130	32.4850	1.3855e+006
	T = 10	8.9220	40.8290	5.4083e+030
	T = 50	15.4380	47.4060	4.8989e+153
	T = 100	19.3600	51.0950	2.4344e+307


Benjamin Mathon
email: benjamin.mathon@gipsa-lab.grenoble-inp.fr

Patrick Bas
email: patrick.bas@gipsa-lab.grenoble-inp.fr

François Cayre
email: francois.cayre@gipsa-lab.grenoble-inp.fr

Benoît Macq
email: benoit.macq@uclouvain.be

Optimization of Natural Watermarking Using Transportation Theory

Keywords: H.4.m [Information Systems Applications]: Miscellaneous; D.2.8 [Software Engineering]: Metrics-complexity measures, performance measures Algorithms, Security, Theory Natural Watermarking, Security, Transportation problem, Hungarian method

This article presents a new approach to use secure spreadspectrum watermarking in the Watermarked Only Attack (WOA) framework (the adversary owns several marked contents and wants to estimate the secret key). Because security of watermarking schemes relies on the distribution of marked contents in a private subspace, we use the transportation theory in order to derive a mapping that will minimize the global distortion of Natural Watermarking (NW) while keeping the same security level. Contrarily to another technique presented in [13], the embedding does not depend on computationally demanding pre-computed mapping. Moreover, we show that this method does not reduce the robustness of the classical NW scheme and enables to decrease the WCR by several dBs regarding former implementations. Tests are made on 2000 Gaussian signals.

INTRODUCTION

The transition from analog to numerical data have permitted a best storage and indexation of multimedia contents. However, the expansion of internet and sharing networks have accelerated the piracy on intellectual property. Watermarking is a mean to solve this problem. It consists in embedding information into a digital content (audio, image or video). This mark can be used to link a content to its owner, to allow or not copies of the content or to trace any running copies on illegal sharing networks. It is obvious that this mark must be imperceptible, nobody shall recognize that the content is watermarked. Another constraint is the robustness of the scheme. In the case of copyright protection, the mark must resist common media processing (for example a JPEG compression on images). The last constraint is the security of the scheme which has received more and more attention in the watermarking community [START_REF] Kalker | Considerations on watermarking security[END_REF][7] [START_REF] Cayre | Watermarking security: Theory and practice[END_REF][3] [START_REF] Bas | Practical security analysis of dirty paper trellis watermarking[END_REF]. Security respects the Kerckhoffs' principle: a key is the only unknown parameter of the embedding scheme to an adversary. Security attacks relate to estimation of a part or all the secret key. According to the degree of estimation, the adversary can delete the mark, make it unreadable, or guess an important part of the secret message.

This article deals with the WOA (Watermarked Only Attack) framework: an adversary has access to several marked (with the same private key) contents. The secret key we use are codewords location in a secret subspace [START_REF] Moulin | Data-hiding codes[END_REF]. To embed a content, its projection in the secret subspace must be moved to a decoding region which corresponds to the right message. Security in watermarking is linked with the conditional distribution of the marked contents given the secret key [START_REF] Cayre | Kerckhoffs-based embedding security classes for woa data hiding[END_REF]. We propose in this article an optimization of a secure spreadspectrum technique: Natural Watermarking (NW) [START_REF] Bas | Natural watermarking: a secure spread spectrum technique for woa[END_REF] from the distortion point of view. This new method uses results of transportation theory: it consists in computing the optimal way to match the distribution of host contents to a distribution of marked contents (given the secret key) by minimizing the global square euclidean distance. Section 2 recalls basics on spread-spectrum watermarking schemes and more particularly on the NW modulation. The link between security in the WOA framework and distributions of marked contents in the secret subspace is also presented. Section 3 details how we use results of transportation theory to minimize the distortion. Finally section 4 exhibits experiments on 2000 Gaussian signals to quantify the performance of this new method.

NATURAL WATERMARKING SPREAD-SPECTRUM

Notations

This section lists the different notations used in this article. Data are written in italic fonts, functions are noted in roman fonts and vectors and matrices are set in bold fonts. Vectors are written in small letters and matrices in capital ones. x(i) is the i-th component of a vector x and (x(0), x(1), x(2), . . . ) is the content of a vector x. P δ and f δ denote respectively the cumulative distribution function (cdf) and the probability density function (pdf) of a probability measure δ. σ 2

x denotes the variance of a signal x and .|. denotes the usual scalar product.

Spread-spectrum watermarking: assumptions and definitions

We consider a host Gaussian signal x ∈ R Nv and a binary message m ∈ F Nc 2 we want to hide into x. The secret key we use for the embedding are Nc carriers ui which are generated thanks to a pseudo random number generator (PRNG) seeded with K ∈ N. These carriers come as Gaussian vectors obtained with the PRNG and are further orthogonalized (using Gram-Schmidt procedure) with unit variance in order to provide a basis of the private subspace i.e. ∀i = j , ui|uj = 0. The marked signal y is obtained by adding the host signal with a watermark signal w and it is constructed as follows:

y = x + w = x + Nc-1 X i=0 s(m(i), x)ui, (1)
where s : F2 × R Nv → R is a modulation. As a convention, we set s(1, x) < 0 and s(0, x) > 0. The W CR (Watermarkto-Content Ratio) measures the distortion embedding power and is defined as:

W CR [dB] = 10 log 10 " σ 2 w σ 2 x « . (2)
Decoding is achieved by the normalized correlation z between each secret carrier and the marked signal:

zy,u i = 1 Nv y|ui = 1 Nv Nv-1 X j=0 y(j)ui(j). (3)
The i-th estimated bit equals 1 if zy,u i < 0 and 0 if zy,u i > 0.

We also denote the normalized correlation vector z for a vector y: zy = `zy,u 0 , . . . , zy,u Nc -1 ´.

We measure the performance of decoding using Bit Error Rate (BER); if m denotes the estimated message we have:

BER = card ({i : m(i) = m(i)}) Nc . (5)
If we use the terminology of [START_REF] Cayre | Kerckhoffs-based embedding security classes for woa data hiding[END_REF], we want to reach the "stegosecurity". It means that the distribution of the host contents and the marked contents is the same. We have DKL(X||Y) = 0, with DKL denotes the Kullback-Leibler divergence.

Natural Watermarking

Natural Watermarking (NW) [START_REF] Bas | Natural watermarking: a secure spread spectrum technique for woa[END_REF] is a spread-spectrum technique which enables to keep (possibly) scaled versions of the original zx,u i distributions. This distribution is supposed to have circular pdf, possibly scaled by a factor η ≥ 1 (in order to set distortion):

sNW (m(i), x) = " η(-1) m(i) sign (zx,u i ) -1 " zx,u i , (6)
with:

W CR = 10 log 10 " Nc(1 + η 2 ) Nv -1 « . (7)
If η = 1, NW belongs to the stego-security class [START_REF] Cayre | Kerckhoffs-based embedding security classes for woa data hiding[END_REF].

Carriers estimation and BSS problem

Our works focus on the WOA (Watermarked Only Attack) framework: the adversary has access to No contents (marked with different messages but with the same private key) and he/she knows the full source code of the watermarking algorithm (following Kerckhoffs' principle). The only parameter he/she does not know is the secret key (i.e. the carriers in spread-spectrum schemes). If we write the spread-spectrum formula considering the No contents, we have the matrix relation:

Y = X + W = X + US, (8)
with:

Y ∈ MN v ,No (R) the watermarked signals, X ∈ MN v ,No (R) the host signals, W ∈ MN v ,No (R) the watermark signals, U ∈ MN v ,Nc (R) the carriers, S ∈ MN c,No (R) the modulations of embedded messages.
Security attacks in spread-spectrum schemes consist in obtaining as much information as possible about U. The problem of disclosing U is known as a Blind Source Separation (BSS) problem. Principal Component Analysis (PCA) allows for an adversary to estimate the Nc-dimensional private subspace spanned by the carriers U (the private key), if the embedding alters the covariance matrix of the contents. PCA involves the calculation of the eigenvalues of the covariance matrix of Y (taken column-wise). PCA aims at finding the optimal linear transformation for keeping the subspace that has the largest variance (message embedding increases the variance of the signal in the directions of the carriers). This technique deals with subspace-security and it allows an adversary to tamper with the hidden messages of the contents. We use the "normalized chordal distance" [6][16] to measure the precision of the estimation of the subspace spanned by the carriers. If Û denotes the estimated carriers (after orthogonalization), the normalized chordal distance between U and Û is defined by:

dc = 1 √ Nc Nc-1 X i=0 sin 2 (θi) ! 1/2 , (9)
where θ0 . . . θN c -1 denote the principal angles between U and Û [START_REF] Pérez-Freire | Spread-spectrum watermarking security[END_REF]. The normalized chordal distance dc equals 0 when the subspaces spanned by the both matrices are the same and equals 1 (its maximum) when the matrices are orthogonal.

Previous works: minimization distortion via a pre-computed mapping

In order to guarantee security, we revisit the embedding as a way to push the host content into the right decoding regions while fitting the appropriate statistical model. We aim at using a model-based embedding which can ensure the level of security we want while minimizing the global distortion.

In [START_REF] Mathon | Distortion optimization of model-based secure embedding schemes for data-hiding[END_REF] we have proposed to use a matching M : X → Y between Nm host correlation vectors X and Nm marked correlation vectors Y (over the secret carriers) with a minimal square euclidean distance on average. This triplet (X , Y, M) is called an Nm-map. For each host signal x we want to watermark, we find the nearest neighbour of zx in X (for example zx i 0 ) and we compute zy = M(zx i 0 ), the correlations of the corresponding marked signal in Y. Fig. 1 shows an example of this process. This optimal matching (from the distortion point of view) was to be found by the O(Nm 3 ) Hungarian algorithm [START_REF] Kuhn | The Hungarian method of solving the assignment problem[END_REF]. The previous technique is more detailed in [START_REF] Mathon | Distortion optimization of model-based secure embedding schemes for data-hiding[END_REF] and was performed successfully for secure watermarking of white Gaussian signals, this optimization was enabled to reduce the W CR distortion by at least 2.6 dB while keeping the same security level. However, this method exhibits the following disadvantages:

• the Nm-map (host and marked correlations and optimal matching) must be stored somewhere, it is very restrictive because it depends on the number of bits Nc the user wants to hide,

• the cubic complexity of the Hungarian method can be a real problem in a practical case when Nm grows up.

We propose in the next section a new method for matching host and marked correlations in an optimal way based on transportation theory. Contrarily to the Hungarian method, we consider the probability measures of the host and the marked correlations for the computation of an optimal matching.

THE TRANSPORTATION THEORY AP-PLIED TO SECURE SPREAD-SPECTRUM WATERMARKING

The Monge-Kantorovich transportation problem

The transportation theory can be considered as an engineering problem, defined by Monge [START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF] in 1781: we have piles of clay X ⊂ R Nc with distribution µ we want to move into holes Y ⊂ R Nc with distribution ν (µ and ν being probability measures). We consider a cost function c:

c : R Nc × R Nc → [0, +∞[ (x, y) → c(x, y), (10)
where c(x, y) = h(xy) represents the cost for moving a unit mass of clay from x to y. The goal is to find a bijection T :

X → Y with ν = T µ = µ • T -1 (T pushes µ forward                    X Y          
Figure 1: Hungarian based watermarking (Nc = 2, Nm = 3): the correlations of the signal we want to watermark is zx, its nearest neighbour is zx 1 . Thanks to the optimal matching (from the square euclidean distance point of view), we find zy = zy 0 , the correlations of the marked signal.

to ν) which minimizes the total cost to move the clay to the holes. Formally we search a minimizer T for:

inf T Z R Nc c (x, T (x)) µ(x)dx | ν = T µ ff . (11)
A solution T of this problem is called an optimal transport map. This problem is called the Monge-Kantorovitch problem because it was first formalized by Monge and principal contributions were done by Kantorovitch [10][11]. For Nc = 1, an optimal transport map T for a convex cost function h is given by:

T = P -1 ν • Pµ, (12)
More generally, for any Nc, we have the following theorem [START_REF] Knott | On the optimal mapping of distributions[END_REF]: Theorem 1. Knott-Smith optimality criterion: if µ (resp. ν) represents the distribution of X (resp. Y), sufficient conditions that T minimizes the Monge-Kantorovitch transportation problem with c(x, y) = xy 2 are that: i) T(X ) has distribution ν, ii) the Jacobian matrix JT(x) of T is symmetric and positive semidefinite.

We have seen that security of spread-spectrum schemes relies on the distribution of the correlations z between the marked signals and the carriers. Moreover, we have the analytic expression of this distribution for NW [START_REF] Bas | Natural watermarking: a secure spread spectrum technique for woa[END_REF], we are now able to use the previous theorem in order to match the host correlations to the marked correlations minimizing distortion constraint.

Application to Natural Watermarking

One property of NW is that the distribution of the correlations over the secret carriers before and after watermarking can be modeled by Gaussian distributions (host signals are Gaussian distributed). From [START_REF] Bas | Natural watermarking: a secure spread spectrum technique for woa[END_REF], we have:

zx,u i ∼ N (0, σ 2 x Nv ) = µ. ( 13
)
We consider we embed the constant message m = (0, . . . , 0) for each host content. We have:

zy,u i ∼ e N (0, η 2 σ 2 x Nv ) = ν, (14)
where, if δ = N (0,

η 2 σ 2 x
Nv ), the distribution ν = e N (0,

η 2 σ 2 x
Nv ) is given by its cdf:

Pν (t) =  0 if t < 0, 2P δ (t) -1 if t ≥ 0. ( 15
)
Then an application of the previous results about transportation theory is relevant, we want to find an optimal transport map T0 which matches the distribution of zx,u i to the distribution of zy,u i . The strategy is to use an optimal transport map for each dimension of the Nc dimensional private subspace thanks to Eq. 12. The computation of the optimal transport map T0 for embedding the message (0, . . . , 0) is given by:

T0 : R Nc → R Nc zx → zy = T0(zx), ( 16
)
where

T0 0 B @ zx(0) . . . zx(Nc -1) 1 C A = 0 B @ P -1 ν • Pµ(zx(0)) . . . P -1 ν • Pµ(zx(Nc -1)) 1 C A . ( 17
)
We have:

Pµ(t) = 1 2 " 1 + erf " t √ Nv σx √ 2 «« , (18)
and

P -1 ν (t) = P -1 δ " t 2 + 1 2 « = ησx √ 2 √ Nv erf -1 " t - 1 2 « . (19)
The constructed function T0 is an optimal transport map because it respects the conditions of Th. 1 (the proof is given in appendix). Thanks to the property of symmetry of spread-spectrum, to embed others message than (0, . . . , 0), sign changes must be made on the coefficients of zx in the indices that have undergone symmetries. Afterwards, inverse symmetry must be performed after watermarking in order to embed the correct message m. Formally, the corresponding optimal transport map Tm for any message m is given by:

Tm(zx) = R -1 (T0(Rzx)), (20) with:
R ∈ MN c,Nc (R), R(i, j) =  0 if i = j, (-1) m(i) if i = j. (21)
We can now define a new spread-spectrum technique based on transportation theory, namely Transportation Natural Watermarking (TNW), the corresponding modulation is given by: sTNW (m(i), x) = Tm(zx,u i )zx,u i .

(

) 22
One can see that we have considerably reduced the complexity of the problem of optimal matching: instead of using O(N 3 m ) algorithm, the matching is directly given by Eq. 20. Fig. 2 shows the modulations sNW and sTNW functions of zx,u 0 (Nc = 1, Nv = 512, η = 2, σ 2

x = 1). Contrarily to NW, the modulation of TNW is not piecewise linear. x = 1). Contrarily to NW, the modulation of TNW is not piecewise linear.

EXPERIMENTS

In this section we generate several Gaussian signals and we watermark them with independent messages by using classical Natural Watermarking (NW), Hungarian Watermarking (HNW, see Sec. 2.5) and Transportation Natural Watermarking (TNW). On the one hand we want to show that TNW does not modify the security and the robustness of the classical NW. On the other hand, one shall keep in mind that we are targeting the same security level as with HNW but with significantly lower algorithmic complexity.

Assessments

We generate No = 2000 N (0, 1) Gaussian signals of Nv = 512 components. For HNW, we use an 10000-map constructed with different standardized Gaussian signals. We hide Nc = 10 bits into each signal. We model robustness attacks by adding a Gaussian noise n to y. Attack strength is assessed by means of the W CNR (Watermarked Contentto-Noise Ratio):

W CNR [dB] = 10 log 10 " σ 2 y σ 2 n « . (23)

Security

Distribution of host and marked signals

Fig. 3 shows the distribution of the correlations of the host, NW and TNW signals over the two secret carriers for η = 1. Obviously, distributions of NW and TNW are identical to that of the host distribution. Therefore, we conclude that TNW does not impair the security of the scheme.

Principal Component Analysis

Principal Component Analysis can be used to estimate the subspace defined by the secret carriers. If an attacker can construct a basis of this subspace, he can make embedded messages unreadable (by nullifying the projections on the estimated subspace). More precisely, he recovers an orthogonal base with a PCA: Û of the subspace defined by the carriers U. This estimation depends on the strength of the embedding [START_REF] Pérez-Freire | Spread-spectrum watermarking security[END_REF]. Fig. 4 shows the normalized chordal distance between the secret and the estimated carriers by PCA functions of the parameter η by using NW and TNW.

As can be seen, the estimation of the carriers subspace for TNW closely follows that of NW.

Robustness

Distortion

Tab. 1 quantify the distortion (by means of the W CR) of NW, HNW and TNW for two values of η = 1, 2). As expected, using an optimal transport map enables to reduce the original NW distortion by 3.76 dB when η = 1. Fig. 6 shows the W CR distortion of HNW and TNW functions of Nc. As can be seen the improvement regarding W CR distortion is better when we use TNW instead of HNW when Nc grows up. The reason is that even if the Hungarian method is optimal to find the matching which minimizes the total cost; HNW computes the embedding by looking for the nearest point registered in the pre-computed optimal mapping; the square euclidean distance between one point and the closest point recorded in the host correlations of the used Nm-map grows with respect to the dimension (Nc).

CONCLUSIONS

In this article we have developed an improvement of natural watermarking: transportation theory has been used in order to derive an optimal embedding formula for embedding Nc bits. This new mapping enables us to minimize the global distortion by keeping the same security and the same robustness as the classical scheme. Moreover this embedding method is better (regarding the distortion) and faster (regarding the complexity) than the Hungarian-based watermarking technique proposed in [START_REF] Mathon | Distortion optimization of model-based secure embedding schemes for data-hiding[END_REF]. As can be seen, the improvement of distortion is better for TNW than HNW when Nc grows up. For HNW, we lost optimality of Hungarian method by selecting nearest neighbours of our host correlations in the Nm-map, η = 1.

ACKNOWLEDGMENTS

Pole IAP-VI fund programme. We acknowledge fruitful discussions with Cedric Villani (Ecole Normale Superieure de Lyon) on optimal transportation theory (particularly about the Knott-Smith criterion applied with Nc 1-dimensional transport maps).

Figure 2 :

 2 Figure 2: Modulations sNW and sTNW functions of zx,u 0 (Nc = 1, Nv = 512, η = 2, σ 2x = 1). Contrarily to NW, the modulation of TNW is not piecewise linear.

Figure 3 :

 3 Figure 3: Normalized correlations of host (a), NW (b) and TNW (c) over the first two secret carriers, No = 2000 observations.

Figure 4 :

 4 Figure 4: Normalized chordal distance between the secret and the estimated carriers by PCA functions of the parameter η for NW and TNW. As can be seen, the estimation of the carriers subspace for TNW closely follows that of NW, Nc = 10.

Fig. 5

 5 Fig. 5 depicts the BER (expectation on the 2000 marked signals) w.r.t. the W CNR. For two values of η = 1, 2 we can see that the robustness against AWGN attack is the same for NW and TNW.

Figure 5 :

 5 Figure 5: BER w.r.t. W CNR [dB] for NW and TNW with η = 1, 2, Nc = 10. η = 1 η = 2 NW -14.09 -10.06 HNW -16.74 -12.29 TNW -17.85 -12.79

Figure 6 :

 6 Figure 6: W CR w.r.t. Nc for HNW and TNW modulation.As can be seen, the improvement of distortion is better for TNW than HNW when Nc grows up. For HNW, we lost optimality of Hungarian method by selecting nearest neighbours of our host correlations in the Nm-map, η = 1.

Table 1 :

 1 W CR distortion for NW, HNW and TNW for η = 1, 2, Nc = 10.

	14.09 -10.06
	HNW -16.74 -12.29
	TNW -17.85 -12.79

APPENDIX

In this appendix, we consider µ = N (0, σ 2

x Nv ) and ν = e N (0,

Optimality of the transport map constructed in Eq. 17.

We want to show that T0 0 B @

x(0) . . .

x(n -1)

is an optimal transport map by applying the Th. 1. The condition i) is verified because a signal x is Gaussian iff each coordinate x(i) is Gaussian (we use the separability of the multidimensional Gaussian distribution). For the condition ii), we introduce the Jacobian matrix of T0:

(25) Show that this matrix is symmetric and positive semidefinite. The symmetry is trivial by construction of the transport map T0. A matrix J is positive semidefinite iff the eigenvalues of J are positive or null. So we must prove that `P-1 ν • Pµ

´ (t) ≥ 0. We have:

fµ(t) is positive (pdf) and P -1 ν (t) is a non-decreasing function so `P-1 ν ´ (t) ≥ 0. Then `P-1 ν • Pµ

´ (t) ≥ 0 and the Jacobian matrix of T0 is positive semidefinite. The conditions i) and ii) of the Th. 1 are verified: T0 is an optimal transport map.