Yousri Kessentini
email: yousri.kessentini@univ-rouen.fr

Thierry Paquet
email: thierry.paquet@univ-rouen.fr

Abdelmajid Benhamadou
email: abdelmajid.benhamadou@isimsf.rnu.tn

Off-Line Handwritten Word Recognition Using Multi-Stream Hidden Markov Models

Keywords: Off-Line handwriting recognition, Hidden Markov Models, Latin script, Arabic script, multistream, information combination

In this paper, we present a multi-stream approach for off-line handwritten word recognition.

The proposed approach combines low level feature streams namely, density based features extracted from 2 different sliding windows with different widths, and contour based features extracted from upper and lower contours. The multi-stream paradigm provides an interesting framework for the integration of multiple sources of information and is compared to the standard combination strategies namely fusion of representations and fusion of decisions. We investigate the extension of 2-stream approach to N streams (N=2,…,4) and analyze the improvement in the recognition performance. The computational cost of this extension is discussed. Significant experiments have been carried out on two publicly available word databases: IFN/ENIT benchmark database (Arabic script) and IRONOFF database (Latin script). The multi-stream framework improves the recognition performance in both cases.

Using 2-stream approach, the best recognition performance is 79.8%, in the case of the Arabic script, on a 2100-word lexicon consisting of 946 Tunisian town/village names. In the case of the Latin script, the proposed approach achieves a recognition rate of 89.8 % using a lexicon of 196 words.

Introduction

Despite the growing use of electronic documents in all the economic activities during the last years, the use of paper documents is still playing an important role. This is because the technologies now offer convenient and cheap means to capture, store, compress and transfer digitized images of paper documents, in a transparent mode. However, automating the processing of the huge amount of these particular documents requires specialized reading systems. While many of such systems are already providing good performance for several applications like OCR, Bank checks Readers, Forms readers, etc. [START_REF] Srihari | Handwritten Address Interpretation: A Task of Many Pattern Recognition Problems[END_REF] [START_REF] Knerr | Hidden Markov Model Based Word Recognition and Its Application to Legal Amount Reading on French Checks[END_REF] [START_REF] Amato | Results from a performance evaluation of handwritten address recognition systems for the United States Postal Service[END_REF] [START_REF] Gorski | A2iA check reader: a family of bank check recognition systems[END_REF], the enhancement of performance is still required so as to cope with a wider range of document reading applications. In today's business world Latin script is mostly used, but with the increasing communication among the different communities worldwide more scripts are getting integrated into information systems. In this context, we propose in this work a unified approach for the recognition of Latin and Arabic scripts. As we want the approach to be script independent, the system must proceed without explicit segmentation of handwriting into graphemes. This is because explicit segmentation methods generally rely on script specific rules to find segmentation points. According to the proposed strategy it is therefore mandatory that the system can operate on low level frame features such as directional or pixel densities. In order to achieve good discriminative power with such low level features, we attempt to combine multiple feature streams. Various combination strategies have been proposed in the literature [START_REF] Günter | Ensembles of classifiers for handwritten word recognition[END_REF]. They can be grouped into two broad categories: feature fusion methods and decision fusion techniques.

The first category commonly known as early integration [START_REF] Okawa | Multi-band speech recognition in noisy environments[END_REF], consists in combining the input feature streams by projecting them into a unique feature space, and subsequently use traditional HMM classifier to model the combined observations in the unique feature space (see Figure . 1).

Figure. 1 : Feature combination approach

In contrast, decision fusion, known as late integration [START_REF] Prevost | Combining modelbased and discriminative classifiers: application to hadwritten character recognition[END_REF], consists in combining the single stream classifier outputs (decisions). Different feature representations obtained from the word image are modelled and decoded separately by individual HMM classifiers. The decoded outputs are then combined to get the final text output (see Figure . 2). [START_REF] Bertolami | Early feature stream integration versus decision level combination in a multiple classifier system for text line recognition[END_REF]) compare these two combination methods in the case of offline handwritten text line recognition and shows that both combination methods improve recognition performances compared to any recognisers built from the individual feature streams. Furthermore, in their case, the early integration approach outperforms the decision level combination.

Figure. 2 : Decision fusion approach

A particular method within the decision fusion framework of sequence models falls into the multi-stream hidden Markov model paradigm (see Figure 3). Such an approach has been particularly studied in the domain of Automatic Speech Recognition (ASR) and presents multiple advantages [START_REF] Wellekens | The use of meta-HMM in multistream HMM training for automatic speech recognition[END_REF] (Bourlard and Dupont, 1997):  It offers a mean to merge different sources of information such as acoustic and visual inputs in some applications such as audio-visual automatic speech recognition [START_REF] Dupont | Audio-Visual Speech Modelling for Continuous Speech Recognition[END_REF] and automatic meeting analysis [START_REF] Mccowan | Automatic Analysis of Multimodal Group Actions in Meetings[END_REF].

 It can combine several kinds of independent features.

 The combination can be adaptive: some sources of information can be weighted, or even rejected if they are not reliable.

 The topology of the HMM can be adapted to each source of information.

 It allows asynchronous modelling of streams. Despite these many possibilities, multi-stream techniques have not been studied and applied for handwriting recognition. Only the work reported in [START_REF] Gauthier | Strategies for combining on-Line and off-line information in an on-Line handwriting recognition system[END_REF] and [START_REF] Artières | A Hidden Markov Models combination framework for handwriting recognition[END_REF] have investigated these techniques for on-Line handwriting recognition. The authors investigate the cooperation of on-line and off-line handwriting word recognition and propose a general framework to combine temporal and spatial representation of the signal.

In [START_REF] Kessentini | A Multi-stream Approach to Off-Line Handwritten Word Recognition[END_REF], we have investigated the use of synchronous multi-stream HMM (without stream asynchrony) for the recognition of Latin script. In this preliminary work, tests were conducted on a private database using two contour feature streams.

In [START_REF] Kessentini | A Multi-Stream HMM-based Approach for Offline Multi-Script Handwritten Word Recognition[END_REF], we extended this preliminary work by enriching the feature streams as described in this paper (see section 4.2.1). Also, asynchronous multi-stream HMM were tested for the first time for the recognition of Latin and Arabic scripts using a limited lexicon (less than 500 words) of the IFN/ENIT database and the IRONOFF database. This paper extends these two contributions in various aspects. A comparison of the multistream approach to the early and late fusion is conducted. This experiment shows significant improvement of recognition performance for both Latin and Arabic scripts. New experiments have been conducted on the IFN/ENIT database considering the whole lexicon of 946 town names and including the use of set "e" so as to compare with other recently reported works.

Finally, this paper explores the extension of the multi-stream paradigm to more than two streams. Experiments conducted on the IFN/ENIT as well as the IRONOFF databases using up to 4 streams show improvement of the recognition performance but with a significant increase of the computational load. This paper is organised as follows. Section 2 reviews the multi-stream formalism and presents the training and decoding stages. Generalization to the N-streams model and some complexity problems are discussed in section 3. In section 4, the overall system organization is presented, the stages of feature extraction and modelling techniques are illustrated.

Experimental results are given in section 5 using the 2-stream approach. They show that the proposed system gives promising results both on Latin and Arabic scripts and compares favourably to other published experimentations on the same data. Additional experiments have been carried out on IRONOFF-cheque and IFN/ENIT databases using the N-stream approach. We discuss the contribution and the computational costs of this extension.

Conclusion and future works are addressed in section 6.

Multi-stream statistical framework

The multi-stream framework provides a convenient formalism to combine several information sources, namely feature streams, using cooperative Markov models (see Figure 4).

Formalism

This problem can be formulated as follows: assume an observation sequence X representing a handwritten word to be recognized, is composed of K input streams X k and assume that the hypothesized model M (e.g. a word model) for a word is the concatenation of J sub-unit models M j (j=1,…,J) (e.g., grapheme or character models). Each sub-unit model M j is composed of K models M j k (possibly with different topologies) attached to each of the K input streams. While the K stream subunit models are assumed independent of each other, they are forced to recombine at some recombination states ( on Figure 4). The recognition problem can be formulated as the one of finding the word model M * that maximizes the posterior probability given a sequence of observations X:

X) | P(M argmax M M *    (1)
where  is the set of all possible word hypotheses. Bayes formula gives:

P(X) P(M) M) | P(X argmax M * 
(2) P(X) being independent of the model M it can therefore be ignored for the computation of M * .

When we assume equal prior probabilities P(M) for all possible word hypotheses, then the recognition problem consists in determining the model M * that maximizes the likelihood

P(X|M).

Using subunit decomposition gives:

) M | P(X M) | P(X J 1 j j j    (3)
Assuming that each stream is independent, each sub-model likelihood

) M | P(X j j
can be computed as a weight sum of the K streams likelihood as depicted in equation (4).

J K k k k j j j k = 1 j= 1 P (X | M ) = w P (X | M )   (4)
This equation can be generalized by letting the stream combination be any function of the stream likelihoods and a set of weighting parameters, W. So the equation ( 4) can be rewritten as,

J k k j j j= 1 logP (X | M ) = f(W , P (X | M ), k   (5)
Most of the approaches use a linear weighted combination function of log likelihood as follows:

1 1 log ( | ) log ( | ) J K k k k j j j j k P X M P X M      (6)
In practice, different combination rules have been proposed for multi-stream systems, including linear (sum, product …), non-linear (MLP) or others (maximum, minimum, median…) rules. More details are presented in [START_REF] Hagen | Robust speech recognition based on multi-stream processing[END_REF].

Linear combination of rules requires the estimation of the stream weights according to their relative reliability. Many weighting strategies are proposed in the literature including fixed weights which have to be trained prior to application, and adaptive weights which are estimated during recognition [START_REF] Hagen | Robust speech recognition based on multi-stream processing[END_REF].

After describing the general multi-stream formalism, we present in the next sections how to train and decode such models.

Multi-stream training

Training the multi-stream HMM consists of two tasks: the first task is the estimation of its HMM stream component parameters (mixture weights, means, variances, and state transition probabilities) and the second task is the estimation of appropriate stream exponents.

Maximum likelihood parameter estimation by means of the Expectation Maximization (EM) algorithm can be used in a straightforward manner to train the first set of parameters. This can be done in two ways: either by training each stream component parameter set separately, based on single-stream observations, and subsequently combine the resulting single-stream HMMs, or train the entire parameter set (excluding the exponents) at once using the bimodal observations. In this work, the first method is used because computationally it is less complex to train the single HMMs instead of the product HMM.

In order to model the Latin characters, we built 26 uppercase character models and 26 lowercase character models). In the case of Arabic characters, we built up to 159 character models. An Arabic character may actually have different shapes according to its position within the word (beginning, middle, end word position). Other models are specified with additional marks such as "shadda". In both Latin and Arabic script, each character model is composed of 4 emitting states. The observation probabilities are modelled with Gaussian Mixtures (3 per state).

Embedded training [START_REF] Rabiner | An Embedded Word Training Procedure for Connected Digit Recognition[END_REF] is used where all character models are trained in parallel using Baum-Welch algorithm applied on word examples. The system builds a word HMM by concatenation of the character HMM corresponding to the word transcription of the training sample.

The final training step concerns the optimization of the stream weights. It is shown in [START_REF] Potamianos | Discriminative training of HMM stream exponents for audio-visual speech recognition[END_REF] that the maximum-likelihood estimation of the weights fails.

Indeed, maximising the likelihood with respect to the weighting factor  yields to the selection of a stream with the highest likelihood (hence

0   or 1   ).
The authors also

show that additional constraints on the weight can yield to a satisfactory solution. Two different methods have been investigated in this work and will be detailed below. The stream combination weights are estimated using two different strategies.

-

Equal combination weights

This strategy consists in attributing equal weights to the various streams. Its main advantage is that no data is needed for estimation of the weights and no extra time for the calculation of the weights has to be expended.

-Relative frequency weights [START_REF] Hagen | Robust speech recognition based on multi-stream processing[END_REF] Employing the forced segmentation given by the multi-stream decoding algorithm, the ratio between the number of times an expert (ie. a stream) performs best for a given character, and the number of times this character occurs in the database is computed.

, kj k j j n n  
where , kj n is the number of training frames for which expert k has the largest probability, over all experts, for character j, and j n is the number of frames for character j in the training data.

In this work we noticed that both of the weighting strategies perform similarly.

Multi-stream decoding

During recognition the best word model M * that maximizes P (X | M ) has to be determined.

Two solutions have been investigated in the literature:

 Recombination at the HMM state level: Although it does not allow for asynchrony or different topologies of the stream models, it is pretty simple to implement and amounts to perform a standard Viterbi decoding [START_REF] Forney | The Viterbi algorithm[END_REF] in which local probabilities are obtained from a linear or nonlinear combination of the local stream emission probabilities.

 Recombination at the sub-unit level: it can force the streams to be synchronous where synchrony is required at the end of character models and can allow for asynchrony where asynchrony allows the sub unit states to be independent of each other. It requires a more sophisticated decoding procedure than the Viterbi search. Two different algorithms have been proposed to solve the problem of decoding in this case:

 Two level dynamic programming [START_REF] Sakoe | Two-level DP matching -A dynamic programming-based pattern matching algorithm for connected word recognition[END_REF]: Here the decoding takes place in two steps. A first dynamic programming process is applied at the sub-unit level and each sub-model is then scored against arbitrary portions of the frame data. Secondly, sub-models are merged together in order to find the best overall score.

 HMM-recombination (Bourlard and Dupont, 1997): It is an adaptation of the HMM decomposition algorithm [START_REF] Varga | Hidden Markov model decomposition of speech and noise[END_REF]. The HMM decomposition algorithm is a time-synchronous Viterbi search that allows the decomposition of a single stream (speech signal) into two independent components (typically speech and noise). In the same spirit, a similar algorithm can be used to combine multiple inputs stream into a single HMM model.

It was shown in [START_REF] Dupont | Etude et développement d'architectures multi-bandes et multi-modales pour la reconnaissance robuste de la parole[END_REF] that both algorithms are equivalent. In this work, we choose to use the HMM-recombination algorithm which is described below. With this choice, a simple pre-processing step that consists in building the product HMM is introduced prior to using a classical Viterbi decoding algorithm. Using the two-level dynamic programming approach would require much more development.

HMM-recombination algorithm

The principle of the multi-stream HMM is to model independently each stream between two pre-determined synchronization points, using multiple single-stream HMMs (two in this study). In this study, the synchronization states are the character boundaries (see Figure 5).

Decoding based on this integration method would require the computation of the best state sequence for both streams verifying the synchronous recombination rule at the same time. To avoid the computation of two best state paths, the model can be formulated as a composite or product HMM (see Figure 6) where each state is built by merging a K-tuple of states from the K stream HMMs (here, K=2). The topology of this composite model is defined so as to represent all possible state paths given the initial HMM topologies. Decoding under such a model requires computing a single best path using the well known Viterbi decoding algorithm. In this model, the observation log-likelihood conditioned on a composite state a-A, that is composed of the stream states a and A respectively, is given using a weighted sum of loglikelihood combination function by:

1 2 1 2 log P(X (t), X (t) | a -A) = ω log P(X (t) | a) + (1 -ω ) log P(X (t) | A)
where 12 X (t) (sim ilarly, X (t)) is the observation vector corresponding to stream 1 (similarly, stream 2) and  the reliability of stream 1 ( 01   ).

The transition probabilities of the product HMM are derived from the transition probabilities of the 2 single stream HMMs assuming independence of the models between 2 recombination states. For example,

1 2 P (a -B | a -A ) = P (a | a)× P (B | A )

Generalization to more than 2 streams

Combining more than two streams is a trivial extension of the approach. The product HMM can easily be computing by combining N single-stream HMMs as follows. Let   i i i j S s ; j 1,...n  be the set of states of the single-stream model i. Then the product HMM is defined in the product state space denoted  S as follows:

    1 N 1 2 N 1 2 N i j k l 2 S S S ... S s s , s ,..., s ; j 1,...n ,k 1,...n ,l 1,...n          
In this model, the observation log-likelihood conditioned on a composite state i S  is computed using a weighted sum of the single-stream log-likelihood by:

         N 1 2 N i i i 1 log P X ( t ), X ( t ),..., X ( t ) | S log P X t | f S          
where f is a function that identifies the state in the  th stream HMM that accounts to product

state i S  as follows:    i k f S s    
The transition probabilities of the product HMM are derived from the transition probabilities of the N single stream HMMs assuming conditional independence of the N stream state sets:

          N i j i j 1 P S | S P f S | f S          
The number of states per character in the composite HMM is equal to the product of the number states per single-stream HMM character. For instance, a three stream character HMMs, each made of 4 states, would result in a 444 = 64 states composite HMM etc…

As we can see, complexity is a major concern when dealing with multiple streams. A straightforward solution can be envisaged by reducing the number of states per character model. This possibility has been investigated in the following experiments (see section 5.4).

Some other possibilities for complexity reduction have been left for further developments and will be discussed in see section 5.4.

Discussion

Now some comments on using the multi-stream framework as opposed to the classical single stream approach can be given. First, the multi-stream framework gives the possibility to combine feature streams asynchronously. This means that some features detected at the same location in the image can be desynchronised during decoding. Such possibility is impossible with standard single-stream HMM and is especially suited when dealing with 2D features as is the case in this study dedicated to Off-Line recognition. Second, because a multi-stream HMM can be implemented with the product HMM of the N stream models, the formalism is applicable to any kind of topology, ranging from left/right models to fully connected ergodic models. In addition, we can notice that the formalism can easily be extended to variable duration state models straightforwardly [START_REF] Rabiner | A tutorial on Hidden Markov Model and selected applications in speech recognition[END_REF]. Indeed, assuming independence of each stream model, the duration probability of the product state i S  is the product of the duration probability of each state that account for i S  as follows:

     N S f S d p d p i i 1 ) ]( [ ) ( ) (  
The product HMM framework also allows the multi-stream model to be adapted to a new dataset in a straightforward way by using standard adaptation approaches. Notice however that adaptation of a product HMM is computationally demanding due to its complexity especially with a large number of streams. This particular aspect of multi-stream adaptation would require more research efforts to receive efficient low complexity implementations.

A script-independent recognition system

Earlier there have been a few efforts behind designing multi-script recognition systems [START_REF] Malaviya | Multi-script handwriting recognition with FOHDEL[END_REF] [START_REF] Chaudhuri | An OCR system to read two Indian language scripts: Bangla and Devnagari (Hindi)[END_REF]. However, these systems do not use HMMs, the framework that have lately been successfully applied to develop many handwriting recognition systems. An HMM based framework offers several advantages allowing us to propose a unified script-independent recognition system. These advantages are mainly due to automatic training of character models on non-segmented words (embedded training), and the segmentation-free recognition paradigm that fits particularly well to a script-independent approach. The proposed system (see Figure 7) is based on a multi-stream HMM.

In the following sub-section, we describe the different stages of our approach. In the first step, pre-processing is applied to the word image. Two types of features are considered in this work: (i) contour based features and (ii) density based features. Contour based features are extracted from the lower and the upper contours, and density based features are computed on two different sliding windows with varying width (see Figure 8). Therefore, each feature type (contour or density feature) defines two feature streams representing the input word image.

Each stream model is then separately trained using Baum Welch algorithm. It can be noted that for assigning weights to the multi-streams, both the equal weight and relative frequency strategies (refer section 2.2) perform similarly. The results reported in section-5 are based on using equal weight strategy. The last step is the recognition during which the HMM models are simultaneously decoded according to the multi-stream formalism presented above.

Pre-processing

Pre-processing is applied to word images in order to eliminate noise and to simplify the procedure of feature extraction. It is worth noticing that these pre-processing methods are script independent.

 Normalization: In an ideal model of handwriting, a word is supposed to be written horizontally and with ascenders and descenders aligned along the vertical direction. In real data, such conditions are rarely respected. We use slant and slope correction so as to normalize the word image [START_REF] Kimura | A Lexicon Directed Algorithm for Recognition of Un-constrained Handwritten Words[END_REF].  Base line detection: Our approach uses the algorithm described in [START_REF] Vinciarelli | Offline Recognition of Unconstrained Handwritten Texts Using HMMs and Statistical Language Models[END_REF] based on the horizontal projection curve that is computed with respect to the horizontal pixel density (see Figure 8). Baseline position is used to extract baseline dependent features that emphasize the presence of descenders and ascenders.

Features extraction

An important task in multi-stream combination is to identify features that carry complementary information. In order to build the feature vector sequence, the image is divided into vertical overlapping windows or frames. The sliding window is shifted along the word image from right to left (in case of Arabic words) or left to right (in case of Latin words) and a feature vector is computed for each frame. features. This kind of features has proved to be discriminative for off-Line handwriting recognition, especially for Latin script [START_REF] Kimura | A Lexicon Directed Algorithm for Recognition of Un-constrained Handwritten Words[END_REF]. The second one is based on foreground (black) pixel densities and has been tested only on Arabic script in [START_REF] El-Hajj | Arabic handwriting recognition using baseline dependent features and Hidden Markov Modeling[END_REF]. In this study, we will thus evaluate the discriminative power of these features using two different scripts (Latin and Arabic scripts).

Contour features

These features are extracted from the word contour representation. Each word image is represented by its lower and upper contours (see Figure 9). A sliding window is shifted along the word image, two parameters characterize a window: window width (8 pixels) and window overlap between two successive positions (5 pixels). For each position of a window, we extract the upper contour points (similarly, the lower contour points). For every point in this window, we determine the corresponding Freeman direction [START_REF] Freeman | On the encoding of arbitrary Geometric Configurations[END_REF] and the directions points are accumulated in the directional histogram (8 features).

Word image contour

Word image upper contour Word image lower contour -Lower contour (see corresponding p points as marked red in Figure 10).

-Interior contour on closure (see blue points in Figure 10).

-Upper contour (see yellow points in Figure 10).

-No point found (see green points in Figure 10).

The black points in Figure 10 represent the lower contour. The histogram of the four kinds of points is computed in each window. This second feature set provides additional information about structure of the contour like the loops, the turning points (for instance, one in letter 't'), the simple lines, and the end points on the word image (altogether, four different features).

Word image Image contour

The third feature set indicates the position of the upper contour (similarly, lower contour) points in the window. For this purpose, we localize the core zone of the word image. More precisely, we extract the lower and upper baselines of word images. These baselines divide the image into 3 zones: 1) a middle zone, 2) the lower zone, 3) the upper zone. This feature set (3 features) provides additional information about the ascending and the descending characters, which are salient characteristics for recognition of the Latin script, as well as of the Arabic script. Hence, in each window we generate a 15-dimensional (8 features from chain code, 4 features from the structure of the contour and 3 features from the position of the contour) contour (for upper or lower contour) based feature vector.

Density features

Here we recall the definition proposed in [START_REF] El-Hajj | Arabic handwriting recognition using baseline dependent features and Hidden Markov Modeling[END_REF]. From each frame 26 features are extracted for window of 8-pixel width (and 32 features for window of 14-pixel width).

There are two types of features: features based on foreground (black) pixel densities, and features based on concavity. In order to compute some of these features (for example, f 2 and f 15 as described next) the window is divided into cells where the cell height is fixed (4 pixels in our experiments) see Figure 11.

Let H be the height of the frame in an image, h be the fixed height of a cell, w the width of a frame (see figure 11). The number of cells in a frame n c is equal to n c =H/h. Let n t (i) the number of foreground pixels in cell i in frame t, and b t (i) the density level of cell i in frame t, then:

b t (i)=0 if n t (i)=0 else b t (i)=1
Let LB be the position of the lower baseline, UB be the position of the upper baseline. For each frame t, the features are the following:

     c n i t t i b i b f 2 2 ) 1 ( ) (
-f 3 : difference in y position of gravity centers of foreground pixels in the current frame and in the previous one.

1 ( ) ( 3    t g t g f
where g is computed as

     H j t H j t j r j r j t g 1 1 ) ( ) ( .
) ( where r t (j) denotes the number of foreground pixels in the j-th row in the frame t, -f 4 to f 11 : densities of black pixels for each vertical column of pixels in each frame (note that the frames here are of 8-pixel width).

-f 12 : vertical position of the center of gravity of the foreground pixels in the whole frame with respect to the lower baseline. The result is then normalized by the height H of the frame.

Stream direction (right to left

) H h w H LB t g f  
) ( 12-f 13 -f 14 : density of foreground pixels over and under the lower baselines for each frame. 15where k is the cell that contains the lower baseline.

     c n k i t t i b i b f ) 1 ( ) (
-f 16 : zone to which the gravity center of black pixels belongs with respect to the upper and lower baselines (above upper baseline, a middle zone, and below lower baseline).

f 17 to f 26 : five concavity features in each frame and another five concavity features in the core zone of a word, that is, the zone bounded by the upper and lower baselines. They are extracted by using a 33 grid as shown in Figure 12. we will see that it can be useful for the recognition of Latin script.

Experiments and Results

The proposed recognition system is script independent. It proceeds without explicit segmentation of handwriting into graphemes. This is because explicit segmentation methods generally rely on script specific rules to find segmentation points. According to the proposed strategy it is therefore mandatory that the system can operate on low level frame features such as directional, contour or pixel densities. In order to achieve good discriminative power of such low level features, the proposed approach is based on the multi-stream paradigm that provides a interesting way of combining individual feature streams.

To evaluate the performance of our recognition system, experiments have been conducted on two publicly available databases: IFN/ENIT benchmark database of Arabic words and IRONOFF database for Latin words (French and English). In all experiments, word-level recognition accuracies have been computed.

IFN/ENIT database

The IFN/ENIT [START_REF] Pechwitz | IFN/ENIT-DataBase for Handwritten Arabic words[END_REF] contains a total of 32,492 handwritten words of 946

Tunisian town/villages names written by 411 different writers. Some town/village names occur in the database with slightly different writing style according to the presence or absence of "shadda" for example. It follows that our lexicon is made of about 2,100 valid entries. Four different sets (a, b, c, d) are predefined in the database for training and one set (e) for testing.

Table1 shows the experimental results of the performance of our recognition system using 4 different single streams (upper contour, lower contour and density with two windows varying in their widths; Density1 and Density2 correspond to the windows of 8-pixel and 14-pixel widths, respectively) as a function of the size of the list of word hypothesis. The best recognition rate is 70.5 % obtained using upper contour feature. From these results it appears that the upper contour is significantly better than the three other feature stream for the recognition of Arabic script. To improve the performance given in Table 1, we try to combine the 4 single streams according to the multi-stream formalism. Six possible pairs of streams can be formed. The recognition results of the 2-streams HMM are presented in Figure 13. In all these experiments, we notice that the multi-stream approach improves the performance obtained with any of the single stream HMM. The best 2-stream recognition rate is 79.6% in Top 1 and is obtained by combining upper contour and density2 features. The gain is 9.1% compared to the best single stream recognition rate. Also combining density and contour feature streams performs better than combining 2 contour streams or 2 density streams.

To compare the multi-stream approach to the standard combination strategies namely fusion of features and fusion of decisions, we report the best 2-stream result obtained by combining the features corresponding to the upper contour and the Density2. As shown in Table 2, the multi-stream approach performs better than the two others standard combination strategies of each individual stream model. Notice that the late decision fusion is performed here using equal weights of the two stream likelihoods using the same combination operator as in the multi-stream approach.

In order to compare our results to the works presented in the literature, we report on table 3 the results obtained on the same database (learning sets: a, b, c, d; test set: e) during the international competition in Arabic handwriting recognition systems at ICDAR 2005 [START_REF] Märgner | Arabic handwriting recognition competition[END_REF]. In this table we have reported the performance of the best proposed multi-stream system which shows a gain of 3.67% for the Top 1 solution compared to the winner system of the ICDAR'05 competition. These results should also be compared to those obtained during the last ICDAR'07 contest [START_REF] Märgner | Arabic Handwriting Recognition Competition[END_REF]. We report in table 4 the results obtained during this competition using a learning sets "a, b, c, d, and e" and sets "f" and "s" for testing. Our system shows the best performance on set "s" and the third best performance on set "f". words.

Table 5 reports the performance of our recognition system using 4 different single streams presented earlier in section 4.1 (upper contour, lower contour and density with two window width). One can notice that the lower contour stream performs significantly poorer than the other three features streams, which perform similarly. This can be explained by the fact that lower contours in Latin words are not generally very informative. Compared to the single stream performance on Arabic script where the upper contour performs significantly better than the other three streams, no single stream feature significantly outperforms the others streams. In order to compare the multi-stream approach with the standard combination strategies, we report in Table 6 the best 2-stream result obtained by combining upper contour and Density1 single streams. We notice that the multi-stream approach performs better than the standard combination approaches. The gain is 3.8 % in Top 1. To compare our results to the related works, we report in table 7 the results obtained on the same database. We notice that our system provides good results on Latin script by achieving the second best performance reported on this database. During these experiments, the multi-stream based system have shown significant recognition results for both Latin and Arabic scripts. The comparison of the multi-stream performance to the classical combination strategies namely, fusion of features and fusion of decisions shows the superiority of the multi-stream approach. We investigate in the next section how the multi-stream approach performs using more than two streams.

Results with N-streams

The previous results have demonstrated the interest of using a multi-stream approach for the recognition of both Arabic and Latin scripts. Following this direction, we are now interested in assessing the superiority of this framework when using more than two streams.

As already stated in section 2.5, extension to N-streams is straightforward but significantly increases the computational cost. In fact, the complexity of the Viterbi algorithm is ). As a consequence, the complexity exponentially increases with the number of streams. This is actually a severe limitation and this is why we have experimented this framework using small lexicons. We have conducted some experiments on IRONOFF-cheque database, a subset of the IRONOFF database, which is made of only French cheque words (30 word lexicons). Table 8 shows the 2-stream, 3-stream and 4-stream recognition results. We notice that adding one feature stream generally improves the performance especially when the single stream HMMs are less accurate. Nevertheless it considerably slows down the decoding process. During these experimentations we have also investigated the influence of reducing the number of states per character on the recognition performance. The use of 3 states per character model (3s/c) speeds up the decoding procedure but slightly decreases the recognition results as shown in table 8. Once again we observe a slight increase of the recognition performance by increasing the number of streams but for a significant increase of the complexity of the decoding process. More investigations should be conducted so as to reduce the complexity of the N-stream decoding process. This could be studied by introducing some of the known techniques that have been already proposed for large lexicon decoding [START_REF] Koerich | Lexicon-driven HMM decoding for large vocabulary handwriting recognition with multiple character models[END_REF]. Various strategies should be tested such as beam search decoding, pruning, A* search etc… together with the use of a lexical tree search. It should be noted that due to complexity reasons, the N-stream experiments have been carried out on limited lexicons. Choosing a particular limited lexicon may biased the results.

For this reason, the random selection has been chosen on the IFN/ENIT database. On the contrary, because there already exists a small bank-cheque lexicon on the IRONOFF database, random election of dictionaries has not been used. This is why the experimental settings on the IRONOFF database are probably biased by this bank-cheque lexicon. This may explain why the best 2-stream results are obtained in this case using a different stream combination (stream 2-3) than those obtained in the previous experiment when using the whole IRONOFF lexicon (stream 1-3). This phenomenon is not observed on the N-stream IFN/ENIT experiment, probably due to the random selection of the lexicon that makes the performance more comparable to those obtained in the previous experiments with the whole IFN/ENIT lexicon.

As a conclusion, it must be stated that the main purpose of these experiments was to establish the contribution of additional streams to the recognition performances. Regarding this investigation, the experimental results show in both cases (IRONOFF and IFN/ENIT) the same moderate improvement of recognition performance when using additional streams.

Conclusion and Perspectives

This paper presents a multi-stream HMM-based approach for off-line handwritten word recognition. The proposed system is script independent. It proceeds without explicit segmentation of handwriting into graphemes and makes use of low level feature sets (directional features and colour density features) irrespective of the scripts. Features are then combined according to the multi-stream paradigm. The developed system has been tested on two publicly available databases: the benchmark database IFN/ENIT (for Arabic script) and IRONOFF database (for Latin script). For both scripts the results show significant improvement while using a multi-stream approach. The comparison of the multi-stream performances to the classical combination strategies namely, fusion of features and fusion of decisions shows the superiority of the multi-stream approach. Moreover, the proposed recognition system provides significant results comparable to the best results reported in the literature on both databases. Future works will consist of testing a new combination rules especially the non-linear ones and a joint training method using a composite HMM. We will be also focused on proposing a unified multi-script recognition system operating independently of the script nature.

Exploring the Transferable Belief Model (TBM) [START_REF] Ph | The Transferable Belief Model[END_REF] to improve the combination rules and test a credal HMM [START_REF] Ramasso | Forward-Backward-Viterbi procedures in the Transferable Belief Model for state sequence analysis using belief functions[END_REF] which combines the generality of TBM and mechanisms of HMM are also one of the interesting futures works.

We are also interested in combining other kinds of features providing more complementary information to further improve the results using a N-streams approach. Feature stream selection should also be addressed using this kind of multi-stream framework.

Figure 3 :

 3 Figure 3: 2-Stream combination approach

Figure 4 :

 4 Figure 4: General form of a k-stream model with recombination states between sub-unit models

Figure 5 :

 5 Figure 5 : Example of a multi-stream HMM with 2 streams and 4 states in each character model.

Figure 6 :

 6 Figure 6 : The corresponding product (composite) HMM.

Figure 7 .

 7 Figure 7. Methodology for the 2-stream training and decoding



  Contour smoothing: Smoothing eliminates small blobs on the contour.

Figure 8 .

 8 Figure 8. Upper and lower baselines detection

Figure 9 .

 9 Figure 9. Word image contours

Figure 10 :

 10 Figure 10: Contour feature extraction

Figure 11 .

 11 Figure 11. Word image divided into vertical frames (here without overlap)

)

  of transitions between two consecutive cells of different density levels above the lower baseline.

Figure 12 .

 12 Figure 12. Five types of concavity configurations for a background pixel P

Figure 13 :

 13 Figure 13: The IFN/ENIT 2-stream recognition performances

Figure 14 :

 14 Figure 14 : The IRONOFF 2-stream performances

O

  (T.(C.N) 2 ). N being the number of states per character model, C is the number of characters in the model and T the number of distinct observation frames per state. In the product HMM, the number of states increase to C.N S with S being the number of streams. Therefore, the complexity of the Viterbi algorithm becomes O(T(C.N S ) 2 )= O(T.C 2 .N 2S

Table 1 .

 1 IFN/ENIT recognition performance using single stream features

	Models	1	2	TOP	5	10
	1) U. contour	70.5 78.6 86.3 90.4
	2) L. contour	63.5 73.1 82.6 86.4
	3) Density1	65.1	73	80.6 83.2
	4) Density2	68.7 78.1 83.3 86.9

Table 2 .

 2 Multi-stream results vs. decision and feature fusion approaches by combining upper contour and

	density2 features
	Models	1	2	TOP	5	10
	2-stream	79.6 85.7 91.6 94.5
	Decision Fusion	75.4 83.2 89.5 92.2
	Feature Fusion	74.1 82.6 88.4 90.8

Table 3 :

 3 Comparison with other word recognition systems which are presented in (Mätgner et al., 2005): recognition results in % with the IFN/ENIT dataset e (6033 images)

	System	Top 1	Top 5	Top 10
	ICRA	65.74	83.95	87.75
	SHOCRAN	35.70	51.62	51.62
	TH-OCR	29.62	43.96	50.14
	UOB	75.93	87.99	90.88
	REAM	15.36	18.52	19.86
	ARAB-IFN	74.69	87.07	89.77
	Proposed system	79.6	91.6	94.5

Table 4 .

 4 Recognition results in % of correct recognized images on references d, e[START_REF] Märgner | Arabic Handwriting Recognition Competition[END_REF].The IRONOFF database[START_REF] Viard-Gaudin | The IRESTE On/Off (IRONOFF) Dual Handwriting Database[END_REF] contains a total of 31,346 isolated word images from a 196-word lexicon (French and English). Although the database contains both on-Line and off-line information of the handwriting signals, only the off-line information is used for our experiments. The offline handwriting signals are sampled with spatial resolution of 300 dots per inch (DPI), with 8 bits per pixel (256 gray levels). For the experiments reported in this paper, we have used a training set of 20,898 words and a test set of 10,448

	(ID 01: MITRE; IDs 02-04: CACI; ID 05: CEDAR; ID 06: MIE; IDs 07-08: SIEMENS; IDs 09-12: UOB-ENST; ID
	13: ICRA; ID 14: PARIS V)
	ID	set f	set s
	01	61.70	49.91
	02	11.95	8.01
	03	15.79	14.24
	04	14.28	10.68
	05	59.01	41.32
	06	83.34	68.40
	07	82.77	68.09
	08	87.22	73.94
	09	79.10	64.97
	10	81.65	69.61
	11	81.93	69.93
	12	81.81	70.57
	13	81.47	72.22
	14	80.18	64.38
	Proposed system	82.09	74.51

Table 5 .

 5 IRONOFF recognition performance using single stream features

	Models	1	2	TOP	5	10
	1) U. contour	81.2 86.8 91.4 94.6
	2) L. contour	69.8 78.4 88.6 93.4
	3) Density1	81.8 87.6 92.6 95.6
	4) Density2	80	86.2 92.8 95.8
	Similar to the experiment conducted on the IFN/ENIT database, we have conducted several
	experiments using the 6 possible 2-stream combinations of HMMs. The results are presented
	in Figure 14.

Table 6 .

 6 Multi-stream results vs. decision and feature fusion approaches by combining upper contour and

		density1 features
	Models	1	2	TOP	5	10
	2-stream	89.8 93.8 96.8 98
	Decision Fusion	86 91.4	95	97.4
	Feature Fusion 85.2 90.6 94.4 97

Table 7 :

 7 IRONOFF: Comparison with other word recognition systems

	Authors	Performances (%) TOP 1 TOP 5
	(Tay et al., 2001) system 1	86.6	94.2
	(Viard-Gaudin et al., 2005) system 1	87.4	95.8
	(Viard-Gaudin et al., 2005) system 2	89.8	97
	(Tay et al., 2001) system 2	96.1	99.1
	Proposed system	89.8	96.8

Table 8 :

 8 IRONOFF-Cheque performances performances on the IFN/ENIT (testing on set d, learning on sets a, b, c) are given in Table9using 2 lexicon sizes and 4 states per character model. We use a random selection of dictionaries for each word test sample.

	Models	4 s/c Top1 Top5	3 s/c Top1 Top5
	2-stream 1)2)	86	95.4	87	97.6
	2-stream 1)3)	92.8	99.2	91	99.4
	2-stream 2)3)	94	99.6	89.4	97.8
	2-stream 3)4)	91	99	86.2	97.2
	3-stream1)2)3)	94	99.6	91.6	99.6
	3-stream1)3)4) 92.6	99.2	90.8	98.8
	3-stream2)3)4) 93.4	99	89.4	98.4
	3-stream1)2)4) 93.6	99.8	92	99
	4-stream	94.2	99.2	92.8	99.2

Table 9 :

 9 IFN/ENIT performances on set d)

			Lexicon size
	Models	30		100
		Top1 Top5 Top1 Top5
	1) U. contour	97.8 99.8 95.4	99
	2) L. contour	97.8 99.8 93.4 99.2
	3) Density1	92	97.4 88.4	94
	4) Density2	92.6	98	88.8 94.4
	2-stream 1)2)	98	99.8 95.6 99.8
	2-stream 1)3)	98.6 99.8 98.4 99.4
	2-stream 1)4)	98.7	100 98.6 99.8
	2-stream 2)3)	98	100 96.8 99.8
	2-stream 3)4)	97.6 99.6 96.4 98.6
	3-stream1)2)3) 98.6 99.8 98.4 99.8
	4-stream	99.2	100 98.8	100


Florian Yger

Maxime Berar

Gilles Gasso

Alain Rakotomamonjy

A supervised strategy for deep kernel machine

. In our proposed procedure, the hidden layers are learnt in a supervised fashion based on kernel partial least squares regression. The main interest resides in a simplified learning scheme as the obtained hidden features are automatically ranked according to their correlation with the target outputs. The approach is illustrated on small scale real world applications and shows compelling evidences.

Introduction

Deep architectures represent a prominent research area in machine learning. However, this approach was for a long time penalized by its inability to derive an efficient learning method leading to a model with good generalization skills. In that sense, shallow architecture such as kernel methods [START_REF] Schölkopf | Learning with kernels: Support vector machines, regularization, optimization, and beyond[END_REF] were preferred. Recent attention devoted to deep architectures stems from the so-called deep learning strategy [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF]. This strategy consists in layer-wise unsupervised pretraining followed by a final supervised tuning [START_REF] Bengio | Greedy layer-wise training of deep networks[END_REF][START_REF] Vincent | Extracting and composing robust features with denoising autoencoders[END_REF]. Doing so allows the hidden layers to learn useful internal representations of data.

Another trend of deep learning consists in pre-training the hidden layers of a deep model in a supervised way. This idea pursued in the 90's was declined in [START_REF] Lengellé | Training mlps layer by layer using an objective function for internal representations[END_REF] where internal features were learnt by minimizing a Fisher discriminant criterion. In the same streamline, Cho and Saul [START_REF] Cho | Kernel methods for deep learning[END_REF] has proposed the multilayer kernel machines (MKM) which represents one attempt to bridge the gap between deep architectures and kernel methods. First, these authors have proposed a new family of kernel functions which mimic the behavior of an infinite size deep architecture and turn a neural network into a shallow structure trainable under the convex optimization framework of SVMs. The so-designed kernel was also applied to train the hidden layers of a deep network.

In the procedure highlighted by Cho and Saul, the hidden layers are greedily trained by first applying a kernel PCA projection followed by features selection based on mutual information with class labels and supervised cross-validation. The obtained model proves competitive with usual deep architectures. However, a major drawback is the complication of features selection stage which involves a not easy setup of many hyper-parameters.

The main contribution of our paper it to alleviate this drawback by learning each hidden layer using kernel partial least squares regression (KPLS) [START_REF] Rosipal | Kernel PLS-SVC for linear and nonlinear classification[END_REF]. The interest of our proposal resides in a simplified learning scheme as the features provided by KPLS are automatically ranked according to their covariance with the target outputs. Therefore, we do not need to resort to any extra crossvalidation. As a side effect, the learnt features are more informative as shown by the empirical evidences in Section 4.

Deep Kernel Machines

Assume a discriminative task with samples set

D = {(x i , y i ) ∈ R q × R p } N i=1 .
Deep kernel machines paradigm aims at learning useful features by stacking successive layers. Each layer is issued from a nonlinear projection of output of the previous layer using kernel method. Given a sample x the output of the -th layer is the size m vector h (x) = g(h -1 (x)) where g is a nonlinear function. Specifically, the j-th value of this vector is expressed as

h j (x) = u j , φ h -1 (x) F (1)
with φ a mapping function which projects any input into the feature space F and u j ∈ F. As usual in kernel methods, u j is expanded over the set of basis functions {φ h -1 (x i ) } N i=1 and reads

u j = N i=1 α j i φ h -1 (x i ) , ∀ j = 1, • • • , m (2)
where coefficients α j i are to be sought at each stage of deep kernel machines learning. From equations (1) and (2), we obtain

h j (x) = N i=1 α j i k h -1 (x), h -1 (x i )
where k stands for the kernel function. Once the last layer is built, the resulting features are used to train any classifier.

Beyond the selection of the depth of the architecture, the issues of deep kernel machines concern how to learn the coefficients vector α j in (2) and to specify the size of each hidden unit. The next section is devoted to the presentation of existing solution and our proposed methodology.

Training of Multilayer Kernel Machines (MKM)

As mentioned previously, the architecture is learned layer by layer. The general framework is given in Algorithm 1.

Existing solution [START_REF] Cho | Kernel methods for deep learning[END_REF] is a compound method that combines an unsupervised learning followed by supervised model selection. In contrary, the approach we promote here is simple, efficient and relies on KPLS, a supervised method. Details about these methods are given below.

KPCA based MKM training

To learn each layer , the solution promoted by Cho and Saul [START_REF] Cho | Kernel methods for deep learning[END_REF] computes N principal directions u j from which the best m informative components are selected. According to [2, Chapter 14], u j is an eigenvector and the corresponding

502

ESANN 2011 proceedings, European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning. Bruges (Belgium), 27-29 April 2011, i6doc.com publ., ISBN 978-2-87419-044-5. Available from http://www.i6doc.com/en/livre/?GCOI=28001100817300.

Algorithm 1 Learning of Multilayer Kernel Machine

[optional] Feature selection on raw data for every layer do Compute the principal directions u j (2) and corresponding h j (x i )

[optional] Rank and select the best principal directions Output the representation h (x i ) for each training sample x i end for Apply any off-the-shelves classifier using features of the last layer coefficients α j are solution of the eigenvalue problem λ j α j = K α j where λ j is the j-th eigenvalue of the centered kernel matrix

K = (I -M) K (I -M) with ( K ) ij = k h -1 (x i ), h -1 (x j
) and M a square matrix with entries all equal to 1/N . Once the N principal components (features) are extracted, only the most interesting are retained. For this sake, the algorithm builds class-conditional and marginal histograms to estimate feature's mutual information with the class labels and re-ranks the features accordingly. Then a kNN classifier is applied on the first m features and the validation error is recorded. This step is repeated for different values of k and m in order to select the best m features.

This procedure deems quite complicated as it involves at each layer first an unsupervised learning and then a trial and error approach to find the best features. The method can be simplified considerably in the sense that the crossvalidation step can be rendered optional, depending on the used kernel projection. Our contribution is derived in that simplification sense and relies on KPLS strategy without validation step. This matter is developed hereafter.

Supervised MKM learning

We propose to get rid of the feature selection and embed it in the projection method. Therefore, the optional steps of algorithm 1 are no longer necessary. This is achieved through kernel partial least squares procedure.

Partial Least Squares method (PLS) computes latent variables based upon the eigenvectors of the empirical covariance matrix between the input and the output variables. Extension to kernel framework were proposed and has been used in the context of multiclass classification (see [START_REF] Rosipal | Kernel PLS-SVC for linear and nonlinear classification[END_REF] for full details).

To train each layer of the MKM in a supervised way, KPLS searches for the jth feature h j (x 1 ) • • • h j (x N ) the most correlated with the class labels. Hence, to obtain u j in Equation 2, one solves the eigenvalue problem λ j α j = YY K j α j (for corresponding vector α j ) with K j obtained by deflating j -1 times the matrix K . In the previous relation, Y ∈ R N ×p stands for the class labels matrix. The deflation process aims at making the matrix K j+1 orthogonal to the subspaces generated by the previous features α j , as the information modelled by the feature α j should not be used in the successive features : K j+1 = H j K j H j where H j = I - procedure with the deflated matrix K j+1 . This scheme considerably simplifies informative features extraction. Indeed, as the features generated by KPLS are obtained iteratively and maximize their covariance with the class-labels, a good guess of the importance ranking of the learnt representations is directly obtained through the eigenvalue λ j . Features number selection is achieved by thresholding these eigenvalues (in practice at a value of 10 -15 ). We coded Y in a one-against-one strategy. Finally, as in [START_REF] Cho | Kernel methods for deep learning[END_REF], we built on top of the last layer a Large Margin Nearest Neighbor classifier [START_REF] Weinberger | Distance metric learning for large margin nearest neighbor classification[END_REF].

K j α j α j K j α j K j 2 α j .

Experiments

In this section, we evaluate our learning approach on three multiclass datasets. The first two datasets were created from the MNIST handwritten digit dataset and were used as a benchmark by Cho and Saul [START_REF] Cho | Kernel methods for deep learning[END_REF] for their MKM. Finally, as a proof of concept, our approach is tested on a texture recognition task.

Although any non-linear kernel could be applied, we follow Cho and Saul and use arc-cosine kernel designed to mimic the behavior of an infinite size neural network. Its expression is given by k

n (x, z) = 1 π x n z n J n (θ) with J n (θ) = (-1) n (sin θ) 2n+1 1 sin θ ∂ ∂θ n π-θ sin θ , the angular dependency θ = cos -1 x,z x z
and n the order of the kernel (see [START_REF] Cho | Kernel methods for deep learning[END_REF] for details). Note that due to ill-conditioning and stability issues, we did not explore order n > 1. Finally, let say that all presented results are averages over 10 runs.

Digit recognition benchmark

Datasets mnist-back-image and mnist-back-random 1 are composed of 28 × 28 grayscale handwritten digits which backgrounds were respectively filled with random image patches or random value pixels. These datasets have the particularity to show the largest error rate gap between shallow and deep architectures. Note that for scalability matters, we used subsets of those two MNIST datasets.

In MKM algorithm, the action performed in the first layer consists in a feature selection. The images being centered, any feature selection method would focus on the central pixels. Therefore, the number of features for the first layer was limited to 300. Even if the subsequent layers were also restricted to 300 features, we observed in practice that less features were selected. We compare the approach based on supervised selection of KPCA features [START_REF] Cho | Kernel methods for deep learning[END_REF] with our KPLS features selection using the same experimental protocol. The results were obtained using validation, training and test sets composed respectively of 500, 1000 and 2000 samples. Figure 1a depicts the error rates obtained for an increasing number of layers by 4 deep kernel approaches on mnist-back-image dataset. For this dataset, KPLS approaches clearly outperform the KPCA based approaches and showed a lower and stable variance. On the other MNIST dataset, results highlighted in figure 1b reveal that none of the methods clearly perform better than the others. Still, the variance of the KPLS approaches keeps on being low. Notice that we have implemented the approach in [START_REF] Cho | Kernel methods for deep learning[END_REF] but we may have missed some heuristics to speed up the algorithm. Hence, due to computational issues, we shall only present results on small parts of the two MNIST datasets. As we were unable to run the KPCA based MKM on the whole MNIST datasets, we used subsets of those. This fact explains the difference in term of classification error rate between our results and those reported in [START_REF] Cho | Kernel methods for deep learning[END_REF]. However, we retrieve the global behavior of deep learning KPCA as exposed [START_REF] Cho | Kernel methods for deep learning[END_REF], even if our experiments are in small scale setting.

Application to texture classification

First of all, it should be stated that the MKM methods are not dedicated to texture recognition. They have been applied out of the box on the data in order to illustrate the gap of performance between the approaches. The Brodatz dataset consists of 112 textures of 640 × 640 pixels (as illustrated in figure 2a), among which we selected 29 images 2 to generate our dataset. As it contains 29 classes, a random classifier should achieve a mean error rate around 96.6%. Inspired by the protocol in [START_REF] Mairal | Supervised dictionary learning[END_REF], we split every texture into two non overlapping parts and randomly extracted 16 × 16 patches. Validation, training and test sets contained respectively 500, 1000 and 2000 samples.

On textures dataset, primary pixel selection could not be applied as for MNIST datasets where informations are spatially localized. Hence, the results for the first layer in figure 2b were obtained on raw textures. Here again, the size of layers > 1 was limited to 300.

KPLS with order 0 kernel and KPCA with order 1 perform the best for this dataset with the trophy for KPLS. However, the impact of the kernel choice is completely different from the previous experiment and it may be explained by the numerical instability we observed as the number of layers increases. Here again, KPLS strategies have the lowest variance in error rate.

Conclusion

In this paper, we proposed a simplified and fast learning scheme for multilayer kernel machines based on KPLS. Experimental comparisons with the original KPCA approach on three datasets validate our method. Moreover the empirical computational cost for KPCA approach is higher (although not presented here). However, some behaviors are still unexplained for both schemes and further studies are still required to fully apprehend the impact of the chosen kernel on the model accuracy. Also, among forthcoming issues to be addressed are the application of other projection methods as for instance kernel CCA or combination of KPLS and KPCA and the large-scale extension of our approach.

Fig. 1 :

 1 Fig. 1: Mean classification error rate at each layer of the MKM .

Fig. 2 :

 2 Fig. 2: Experimental results on Brodatz dataset.

  The next feature is extracted by repeating the 503 ESANN 2011 proceedings, European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning. Bruges (Belgium), 27-29 April 2011, i6doc.com publ., ISBN 978-2-87419-044-5. Available from http://www.i6doc.com/en/livre/?GCOI=28001100817300.

ESANN 2011 proceedings, European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning. Bruges (Belgium), 27-29 April 2011, i6doc.com publ., ISBN 978-2-87419-044-5. Available from http://www.i6doc.com/en/livre/?GCOI=28001100817300.

* This work is funded in part by the FP7-ICT Programme of the European Community, under the PASCAL2 Network of Excellence, ICT-216886 and by the French ANR Project ASAP ANR-09-EMER-001.


Cornelia Vacar

Jean-Franc ¸ois Giovannelli

Yannick Berthoumieu

LANGEVIN AND HESSIAN WITH FISHER APPROXIMATION STOCHASTIC SAMPLING FOR PARAMETER ESTIMATION OF STRUCTURED COVARIANCE

Keywords: Stochastic sampling, Monte Carlo Markov Chains, Hessian, Fisher, texture

We have studied two efficient sampling methods, Langevin and Hessian adapted Metropolis Hastings (MH), applied to a parameter estimation problem of the mathematical model (Lorentzian, Laplacian, Gaussian) that describes the Power Spectral Density (PSD) of a texture. The novelty brought by this paper consists in the exploration of textured images modeled by centered, stationary Gaussian fields using directional stochastic sampling methods. Our main contribution is the study of the behavior of the previously mentioned two samplers and the improvement of the Hessian MH method by using the Fisher information matrix instead of the Hessian to increase the stability of the algorithm and the computational speed.

The directional methods yield superior performances as compared to the more popular Independent and standard Random Walk MH for the PSD described by the three models, but can easily be adapted to any target law respecting the differentiability constraint. The Fisher MH produces the best results as it combines the advantages of the Hessian, i.e., approaches the most probable regions of the target in a single iteration, and of the Langevin MH, as it requires only first order derivative computations.

INTRODUCTION

The sampling methods prove to be important tools in signal and image processing, as they enable the exploration of intricate target laws and permit their characterization by computing their statistical descriptors based on the drawn samples.

The performances of these algorithms directly influence the performances of the algorithms that embed them and this is why the speed and the accuracy with which the resulting samples describe the target are very important. We have investigated the performances of two stochastic samplers, Langevin and Hessian adapted MH, in the context of a parameter estimation problem for textured images.

First order derivative-based samplers have also been employed in [START_REF] Choo | People Tracking Using Hybrid Monte Carlo Filtering[END_REF] and in the recent work [START_REF] Kayabol | Adaptive Langevin Sampler for Separation of t-Distribution Modelled Astrophysical Maps[END_REF], but these approaches explore the Hamiltonian MH algorithm, whilst we have studied the Langevin MH. However, our contribution is more significant in the case of the second order derivative-based methods, where, in the present estimation problem, we have provided an improvement, by modifying the Hessian algorithm presented in [START_REF] Qi | Hessian-based Markov Chain Monte-Carlo Algorithms[END_REF] in the sense that the concerns regarding the positiveness of the Hessian matrix are eliminated by its replacement with the Fisher information matrix. The results obtained for the proposed problem are encouraging, confirming the stable behavior and speed performances of the second order derivative-based samplers.

BAYESIAN SETTING

Although the sampling algorithms are able to explore any type of target, not necessarily a probability density, we will establish from the beginning a Bayesian context of a parameter estimation problem. In this setting, we will use the following notations throughout the paper: x is the set of observed data, θ is the set of parameters to be estimated, f (x|θ) is the likelihood, describing the probability of a certain observed data configuration given the parameters, and π represents the target we intend to sample.

An MH algorithm, as described by [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF], [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF], relies on a transition kernel consisting of two ingredients: an arbitrary transition kernel, q(θ, φ), and an acceptance probability:

α(θ, φ) = min 1, π(φ)q(φ, θ) π(θ)q(θ, φ) (1)
where θ is the current and φ is the proposed value for the parameters. The process thus consists in formulating a proposal according to the transition kernel, evaluating the acceptance probability of the proposal and deciding to accept or reject it based on the computed probability.

The functioning of an MH algorithm is:

1. Initialize the iteration counter j = 1 and set θ (0) . 2. Propose a new value φ for the parameter, generated from the density q(θ (j-1) , •). 3. Evaluate the acceptance probability α(θ (j-1) , φ) given by Eq.( 1) and according to this value accept the proposal and update the parameter θ (j) = φ, or reject it and keep the old value for the parameter θ (j) = θ (j-1) . 4. Update the counter and return to step 2 until convergence.

All analyzed samplers are convergent and accurately explore the target, the differences between them being given by the different formulation of the transition kernel and, implicitly, by the expression of the acceptance probability. These differences translate into different evolutions of the samples when exploring the target (and how quickly the samples begin to be drawn from the most representative regions) and in differences between the time needed by each sampler to produce a sample.

Independent MH has a proposition law that does not depend on the current position of the chain, Random Walk (RW) MH proposes an isotropic displacement around the current value, while more advanced methods include a directional component for the proposal, meaning that the proposal is built as follows: add to the current value the directional component and then add the stochastic component, i.e., make an isotropic move around this new value. In the class of directional MH methods lie the Hamiltonian, Langevin and Hessian adapted MH, the last two algorithms being those that were studied here, with an emphasis on the Hessian method, which has been improved by building a new transition kernel based on the Fisher information matrix.

In the following we will present these algorithms, adapted to our parameter estimation problem. In a Bayesian parameter estimation algorithm the a posteriori law for the parameters, given the observations, is the law that dictates the probability distribution law in the parameter space. Thus, this is the law we should explore, i.e., the target, and extract the statistics from the yielding samples. According to Bayes' law, the a posteriori distribution is proportional to the product between the a priori distribution and the likelihood. With the pre-requisite that the a priori will be embedded in the proposition law, the likelihood accurately approximates the posterior distribution law, and thus, we will sample, the Co-Log-Posterior(CLP ), where CLP = -log f (x|θ), will be evaluated as being the target in order to achieve the convexity of the target with respect to the parameters.

SAMPLING ALGORITHMS

The principle of stochastic methods is the sampling of a distribution by considering it as being the limiting distribution of a Markov chain and simulating this chain until equilibrium is reached, as stated in [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. Furthermore, the class of MCMC methods has been narrowed to the Metropolis Hastings (MH) family, due to its adaptation to sampling intricate target laws.

Langevin Adapted Metropolis Hastings

Langevin algorithms are derived from diffusion approximations and rely on the principle of using the information concerning the target density, in the format ∇ log π, in order to build a proposal distribution well-adapted to the problem in question [START_REF] Gelman | Bayesian Statistics[END_REF]. The Langevinbased MH algorithm proposes a random walk-like transition of the form [START_REF] Besag | Bayesian computation and stochastic systems (with discussion)[END_REF]:

θp = θc + ε 2 2 g(θc) + ε N (0, I) (2)
where g(θc) = ∇CLP (θ)|θ=θ c . The acceptance probability can be obtained from (1), for:

q(θp, θc) = exp - 1 2ε 2 θc -θp - ε 2 2 g(θp) 2 (3)
As compared to the non-directional MH methods, the complexity of the Langevin MH is increased due to the form of the acceptance probability and to the necessity of evaluating the gradient for every new proposal. However, the increased computation time per iteration is compensated by the much smaller number of iterations needed to reach convergence.

In regions far from the maximum of probability, the gradient is large (the directional component is dominant), thus the algorithm approaches the high probability regions with high amplitude jumps, while, when near the maximum of probability, the gradient is small, thus the stochastic component is dominant and the region of high probability is explored due to the stochastic component.

Hessian Adapted Metropolis Hastings

This section is dedicated to a sampling method that has seldom been explored and whose presence in the literature is scarce. The directional component of the proposal is in this case formulated using Newton's direction, which, for a quadratic law, indicates the maximum. In [START_REF] Qi | Hessian-based Markov Chain Monte-Carlo Algorithms[END_REF] a version of this sampler has been tested and compared to methods such as Gibbs and optimal marginal data augmentation (DA) samplers on a probit regression problem, proving that the performances of this sampler are superior. The transition for the Hessian sampler is of the form:

θp = θc + ε Σ(θc) g(θc) + N (0, Σ(θc)) (4)
where Σ(θc) = -H(θ) -1 |θ=θ c . The acceptance probability will be obtained from (1), for:

q(θp, θc) = N (θp -θc + ε Σ(θc) g(θc), Σ(θc)) (5)
The advantage of the method is that, for quadratic and quasiquadratic distributions, the regions of high probability are approached in a very small number of iterations (ideally, a single one) and then explored with the contribution of the stochastic part of the proposition, a quadratic law of variance Σ(θc), which is an accurate approximation of the target. However, it is clear that this method is also rather complex, as each iteration translates in the computation of the gradient and the Hessian matrix and the evaluation of the acceptance probability.

Fisher Adapted Metropolis Hastings

Another concern is the need to perform the inversion of the Hessian and in order to avoid matrix inversion problems that may occur if the Hessian is not positive definite, the idea of replacing the Hessian by the Fisher information matrix has been explored in the present work. The Fisher matrix quantifies the amount of information that the observations contain regarding the parameter θ. We have been able to apply such an approximation, as in the case of our problem, we have a great amount of independent observations, thus, a scenario close to the asymptotic case. In such a situation, the a posteriori distribution is consistent and normal of inverse variance equal to the Fisher matrix [START_REF] Philippou | Asymptotic Distribution of the Likelihood Function in the Independent not Identically Distributed Case[END_REF]. Thus, instead of computing the Hessian, the quantities:

Ipq(θ) = Ex|θ ∂ 2 ∂θp∂θq CLP (θ) θ (6)
are used in order to formulate the proposal and to evaluate the acceptance probability. This approximation reduces the complexity of computations as, instead of computing the Hessian, when applying the expectation with respect to the data given the parameters, the term containing second order derivatives becomes null. This fact is due to the particular problem we are addressing, i.e., the model chosen for our textures and to the fact that the coefficients in the Fourier domain have the specific distribution detailed in Section 4. The performances of the algorithm are enhanced by this approximation, as the time per iteration is reduced because no second order derivatives of the target must be computed and all is reduced to first order derivatives computations. The major advantage brought by this innovation is that matrix inversions are always possible and the algorithm is stable, no matter if there are particular cases when the Hessian is not positivedefinite (case that is mentioned, but not dealt with in [START_REF] Qi | Hessian-based Markov Chain Monte-Carlo Algorithms[END_REF]).

TEXTURE CONTEXT

This section presents the application of the previously presented sampling methods to a parameter estimation problem. The principle is that, starting from texture images such as those in Fig. 2, we estimate the parameters that characterize the features of the texture and x ν 0 y parameter space for the four samplers. To the left, observe the sparseness of the accepted samples for Independent MH, for RW MH the evolution step is very small and undirected, while for Langevin MH the proposal is influenced by the gradient. For Hessian (Fisher) MH, the strong probability regions are approached in a single iteration and then thoroughly sampled, as these are the regions most representative for the target. thus, reduce its description to a small number of numerical quantities. This approach is possible due to the particular nature of these textures, i.e., to the fact that they can be modeled by zero-mean, stationary Gaussian Random Fields (GRF) of covariance matrix R [START_REF] Haindl | Parameter Estimation in Gaussian Markov Random Fields[END_REF].

As seen in Fig. 2, this method of modeling the textures can describe a broad class of images, i.e., by the simple variation of a relatively reduced set of parameters, various patterns and characteristics appear.

Considering x as being the N 2 pixel vector of such an image, the model is completely described by its second order statistics and a distribution law of the form:

f (x|R) = (2π) -N 2 /2 det(R) -1/2 exp - 1 2 x t R -1 x (7)
Working under the stationarity hypothesis, the R matrix has a Toeplitz form and can be approximated by a circulant form R = F † ∆F, where ∆ is a diagonal matrix containing the eigenvalues of R and F denotes the Fourier basis matrix of size N × N . After further computation, Eq.( 7) can be rewritten as:

f (x|R) ∝ N n,m=1 s -1/2 nm exp - 1 2 N n,m=1 | • xnm| 2 /snm (8)
In Eq.( 8), |

• xnm| 2 represent the squared moduli of the image's Fourier coefficients and snm are the discretized elements of the PSD. The principle of the model does not consist in a discretization of the PSD, which could lead to aliasing problems, but in the use of the PSD values at discrete positions in order to describe the law Fig. 2: Texture realizations using GRFs with structured PSDs.

governing the elements of the image in the Fourier domain. The previously mentioned form is achieved using a Whittle approximation that enables the evaluation of the CLP in O(N log N ). In fact, the snm elements dictate the distributions of the image's Fourier transform coefficients in the following manner:

• xnm ∼ N (0, snm) (9)
This means that each Fourier coefficient

• xnm has a zero-mean, normal distribution, of variance snm, thus we are in a scenario of independent, but not identically distributed observations.

The PSD of the texture can be described by a wide variety of positive-valued mathematical laws and the snm quantities that intervene in Eq.( 8) follow these laws. Our study is focused on the Lorentzian, Laplacian and Gaussian functions and in the following the Lorentzian form will be used to exemplify the computations:

SL (νx, νy, θ) = γ 1 + (νx -ν 0 x ) 2 w 2 x + (νy -ν 0 y ) 2 w 2 y ( 10
)
where θ = [γ, ν 0 x , ν 0 y , wx, wy] ∈ R 5 is the parameter set that drives the PSD and SL (νx, νy, θ) denotes the continuous PSD at position (νx, νy) for the current value θ of the parameter set. Furthermore, the previously mentioned snm elements are functions of θ and represent the value of the PSD at discrete positions in the frequency domain, as shown by the expression:

snm(θ) = SL (n∆νx, m∆νy, θ) (11)
Eq.( 11) illustrates the fact that the snm elements contain the contribution of the parameter set θ through the mathematical model for the PSD. Furthermore, as shown by Eq.( 9), these elements of the PSD govern the distribution of the texture's Fourier transform. It is thus obvious that the law of dependency of

• xnm on θ is cumbersome, for which reason stochastic methods were employed.

As mentioned in the previous section, the goal of the sampling algorithms is to properly explore the a posteriori law. The directional character of the analyzed samplers would be fully exploited if the CLP had specific convexity properties, and this is the reason for which the change of parametrization λnm(θ) = snm(θ) -1 is proposed in order to render the target law convex with respect to the elements of the PSD. Consequently, the expression of the CLP , convex with respect to λnm(θ), is: The algorithms presented in Section 2 employ the first and the second order derivatives of the CLP , described by:

CLP (θ) = 1 2 N n,m=1 log 1 λnm(θ) + | • xnm| 2 λnm(θ) (12)
∂CLP ∂θp = N n,m=1 | • xnm| 2 - 1 λnm(θ) ∂λnm(θ) ∂θp (13)
∂ 2 CLP ∂θp∂θq = N n,m=1 1 λ 2 nm (θ) ∂λnm(θ) ∂θp ∂λnm(θ) ∂θq + N n,m=1 | • xnm| 2 - 1 λnm(θ) ∂ 2 λnm(θ) ∂θp∂θq (14)

RESULTS

The sampling methods presented in this paper are all convergent and the use of one or another influences only the speed of convergence and the mixing properties of the algorithm, not the actual equilibrium state.

As the theory announced, all methods yield posterior histograms showing that the samplers explore the law and produce samples with the same posterior mean and variance. Consequently, in order to compare the methods, the convergence time and the computational cost is analyzed for each sampler. The improvement brought by the directional methods consists in the fact that they approach the region of important probability in a small number of steps, proposing transitions in the direction of probability increase and all at the same time they thoroughly explore the regions most representative for the target.

In Fig. 3 it can be seen that, even if it requires the repeated evaluation of gradients, the Fisher RW has speed performances superior to RW MH and better than those of Langevin RW. The parameter ε that intervenes in all the tested methods was tuned in such a manner as to obtain an acceptance rate of approximately 24%, which in [START_REF] Gelman | Bayesian Statistics[END_REF] is proven as being optimal in the sense that it provides the best compromise between the amplitude of the transitions (that ensures the samples are not too correlated) and the acceptance rate.

The advantage of the Fisher MH is given by the fact that it approaches the regions of high probability in a single iteration, requiring for this only the evaluation of gradients, as opposed to gradients and hessian matrices for Hessian MH. This approximation is possible due to our texture context that positions us in an asymptotic scenario of independent, but not identically distributed observations. This simplification brought by the method is given by the new form of the proposal, as the second term in Eq.( 14) becomes null, due to the distribution of the Fourier coefficients in Eq.( 9).

CONCLUSION AND PERSPECTIVES

The Fisher sampling method proved to be the best suited for the problem of parameter estimation in the case of the texture with a structured PSD, for several complex models. The presented directional sampling methods are also just as efficient in the case of any distribution law, with the condition for it to be twice differentiable.

One of the most interesting directions of pursuit is to further exploit the Fisher approximation and perform a more detailed analysis regarding the gain that such an approximation brings in stochastic sampling applications for other types of targets.

Fig. 1 :

 1 Fig. 1: Samples evolution in the ν 0x ν 0 y parameter space for the four samplers. To the left, observe the sparseness of the accepted samples for Independent MH, for RW MH the evolution step is very small and undirected, while for Langevin MH the proposal is influenced by the gradient. For Hessian (Fisher) MH, the strong probability regions are approached in a single iteration and then thoroughly sampled, as these are the regions most representative for the target.

Fig. 3 :

 3 Fig.3: The stabilization of the recursive means for several chains in the case of each of the four studied methods. As expected, all the algorithms converge to the same value, but the Hessian (Fisher) MH reaches equilibrium the fastest. Although superior in terms of computation speed per iteration, the Independent MH and Random Walk MH require a longer interval to converge.

ACKNOWLEDGEMENTS

We would like to thank Olivier Féron and Franc ¸ois Orieux, with whom we have discussed some of the aspects presented in this paper.