P Basset

A Amelot

J Vaissière

B Roubeau

Nasal airflow in French Spontaneous Speech

The goal of this paper is to compare the degree of anticipatory and carryover phenomena in the behaviour of the velum during the production of French spontaneous speech, and to compare the data with the same speech material, but read. Airflow through the nose and the mouth were taken as an indirect indication of the velum behaviour. French is a language of special interest because there is a phonological contrast between oral and nasal vowels. The results show a tendency for nasal airflow to start before the nasal and a strong propensity to spread after the nasal. No differences could be established between speaking styles (spontaneous vs read speech) regarding nasal airflow anticipation and carryover.

Introduction

Our work deals with aerodynamic data. Nasal airflow measurements were collected as an indirect indication of the degree of opening of the velopharyngeal port as there is a valid correlation between air flow and opening when the cross sectional area is less than 0.4 cm2 [START_REF] Warren | Effects of size of the nasal airway on nasal airflow rate[END_REF]. The presence of airflow indicates that the velopharyngeal port is open but its absence does not signify that the velopharyngeal port is closed. Furthermore, lowering of the velum starts before velopharyngeal port opening [START_REF] Benguerel | Velar co articulation in French: an fiberscopic study[END_REF]. While aerodynamic measurements are only an indirect way of collecting data on velic behaviour, they have however the advantage of being non-invasive which is essential for studying spontaneous speech. French is of particular interest because there is a phonological contrast between oral and nasal vowels (achieved by velopharyngeal port opening). All former aerodynamic studies for French dealt with read speech [START_REF] Durand | De la formation des voyelles nasales[END_REF][START_REF] Benguerel | Nasal Airflow Patterns and velar Coarticulation in French[END_REF][START_REF] Cohn | Phonetic and Phonological Rules of Nasalization[END_REF][START_REF] Ferguson | Nasal Epenthesis in Hindi[END_REF]. This paper provides new data on spontaneous speech. Because of the phonological contrast, it is expected that velic coarticulation should be of a lesser extent in French (as compared to English, for example). However, as it takes some time for the velum to lower and to rise again (50 ms according to [START_REF] Ohala | Reprinted from Nasalfest Paper from a symposium on nasals and nasalization[END_REF], a minimum amount of coarticulation is expected. It was hypothesised that less coarticulatory phenomena (anticipation and carryover) would happen in read, better articulated, speech, as compared to spontaneous speech. Spontaneous conversations were recorded while the speaker was chatting with the experimenter between periods of recording read speech during the first session. As it takes some time to set up and calibrate the instruments, the speaker was asked to keep the two nasal tubes, the oral mask and the pair of plate electrodes on his (her) neck on during the complete recording sessions. The spontaneous data obtained in this way sound fairly natural.

Experimental procedure

Instrumentation

Nasal airflow was measured with a pressure transducer attached to two tubes inserted in the nostrils and secured by rubber plugs. Oral airflow was measured with a tube inserted in a soft rubber mask that interfered very little with jaw lowering and not at all with lip Basset P., Amelot A., Vaissière J., Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100.

movements. Sound was recorded with a microphone external to the mask. Electroglottographic (EGG) data were also recorded. Eva work station and Phonedit software (Teston and Galinder, 1990) were used.

In order to verify the adequacy of the response time, airflow data were acquired from two different stations (EVA2 and Pcquirer), with the same speaker saying the sentence: "ton tonton tond ton tonton" [tO) tO) tO) tO) tO) tO) tO) ] 'your uncle shaves your uncle'. The synchronisation between the burst on the waveform and a sharp rise in the airflow at the release of the stop was taken as an indication of a good response time of the machine [START_REF] Barry | Co-articulatory airflow characteristics of intervocalic voiceless plosives[END_REF]. There were no major differences between the two devices (see figures 1 and 2). As reported by [START_REF] Benguerel | Nasal Airflow Patterns and velar Coarticulation in French[END_REF], a short phase of negative airflow before a nasal was observed (figure 1). It corresponds to an increase of volume in the nasal cavity at the beginning of the lowering of the velum. It was considered as an indication of adequacy of the instrument. An accurate determination of the response time of the nasal airflow is however a complicated problem. It remains as a task for the future (Ohala, personal communication).

Subjects and Speech material

The total corpus analysed consists of 11 minutes of spontaneous speech uttered by three females (1mn 46s; 1 mn 20 s; 3 mn 18 s) and one male (4 mn 36 s), speaking Standard French. Spontaneous conversations were recorded during breaks between periods of recording read speech. The topics of discussion range from holidays, to previous studies Basset P., Amelot A., Vaissière J., Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100.

on nasalization, and troubles with computers. The conversations were orthographically transcribed. During a second recording session, each speaker was asked to read and record a great part of what he (she) had said spontaneously during his (her) first session. It was then possible to compare the read and spontaneous styles on the same material.

Audio, aerodynamic and glottographic signals were segmented and phonetically transcribed with the aid of spectrographic representation. Priority to the auditory impressions was given in case of conflicting evidence.

Nasal phones at the beginning and at the end of utterances were discarded as respiration can contaminate the data [START_REF] Van Hattum | Air flow rates in normal speakers[END_REF]. Occurrences containing more than two nasals were also discarded when it was impossible to separate anticipatory from carryover influence, as for example in [sE) sA) sE) kA) t fÂA) ] "cinq cent cinquante francs" 'five hundred and fifty francs'.

As shown in figure 3, anticipation means that NAF (Nasal Air Flow) starts before the onset of the nasal, delay that NAF starts after the onset of the nasal and carryover indicates that NAF spreads over the following phone. The following measurements were made (table1):

-The mean duration of anticipation (mean a ) delay (mean d ) and carryover (mean c ).

-The ratio between the mean duration of anticipation and the mean duration of the proceding phone (Dn-1). If the ratio is greater than 1, it means that anticipatory velic opening starts before the onset of the preceding phone.

-The ratio between the duration of delay and the total duration of the nasal.

-The ratio between the mean duration of carryover and the mean duration of the following phone (Dn+1). If the ratio is greater than 1, it means that nasal flow spreads after the end of the phone following the nasal.

Results and discussion

Anticipation

Table 2 presents the distribution of presence and absence of anticipatory nasal airflow across the contexts of a vowel or voiced/voiceless consonant before a nasal vowel, and of an oral vowel before a nasal consonant. Table 3 gives an overview of the total number of items in the four anticipation contexts as well as the relative frequencies of anticipation, synchronous onset and delay of nasal airflow in these contexts.

In V+Vcontext, there is a tendency (18/32, 56%) to anticipate nasalization during an oral vowel. In most cases, NAF is starting during the preceding phone (ratio=1.3). Some perception tests should be conducted to find out if the nasalized preceding phone is still perceived as oral. A nasalised vowel is not identified by French listeners as phonemically nasal because nasal vowels in French differ from oral, not only by the presence of nasal airflow but also by a special articulatory setting [START_REF] Zerling | Phénomènes de nasalité et de nasalization vocaliques: Etude cinéradiographique pour deux locuteurs[END_REF].

Because low vowels are realised with a lower velum [START_REF] Passavant | Ueber die Verschliessung des Schlundes beim Sprechen[END_REF] and because they need a lower velum to be perceived as nasal [START_REF] Benguerel | Perception of vowel Nasalization in French[END_REF][START_REF] Maeda | Acoustic of vowel nasalization and articulatory shifts in French nasal vowels[END_REF], the results for high and low vowels were presented separately. The higher oral pressure during high vowels favours nasal leaks. There was not enough data in the present study to establish a significant difference between the two vowel types.

In C+V there is a clear difference depending on the voicing of C. This is in agreement with the results of [START_REF] Ferguson | Nasal Epenthesis in Hindi[END_REF]. When C is voiced, anticipation prevails (78/100, 78%). In some cases like for example in "C'est bon ?" [sE bO) ] 'is it all right ?' (see figure 6 (i)), the labial voiced stop /b/ preceding the nasal vowel /O) / seems completely nasalized with no burst visible on the spectrogram and with the presence of formants. Out of context, it sounds like a nasal.

When C is unvoiced, anticipation only happens in one third of the cases (31/90, 34%) and there are many cases of delay (in read speech, for stops: 15/32, 46%). In the cases of anticipation, like in the example on figure 6 (ii) "vacances" [vakA) s] 'holidays', nasal airflow seems not sufficient to prevent the building up of oral pressure and the burst is still visible on the spectrogram.

In V+N a large variability was observed. There were as many cases of anticipation as synchronous onset of NAF (75/162, 46%). Table 4 shows the results with both styles merged together. No clear tendency towards anticipation (53%) or synchronous onset of nasal airflow (47%) can be seen. Table 5: Comparison between the two styles, spontaneous and read.

Carryover

Table 6 illustrates the presence or absence of carryover airflow across the contexts of a nasal vowel before an oral vowel or a voiced/voiceless consonant, and of a nasal consonant before an oral vowel. Table 7 gives an overview of the total number of items in the four carryover contexts as well as the relative frequencies of presence and absence of carryover nasal airflow in these contexts.

Basset P., Amelot A., Vaissière J., Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100. While in [-nasal] [+nasal] contexts, no clear difference appears between anticipation and synchronous onset of NAF, or between styles, the results are very different for carryover, even if in V+V context there is only a slight tendency for nasal airflow to spread on the following phone (13/22, 59%). In V+C, in almost all cases, there is a carryover of nasal airflow on the following phone (64/70, 91%, when C is voiced and 159/14, 92%, when C is unvoiced). However, there is a difference in the mean duration and ratio of airflow depending on the voicing. When C is unvoiced, the ratio of the carryover does not exceed 0.8 (the next phoneme is not completely nasalized); the ratio is 1 when C is voiced (the next phoneme is completely nasalized). As shown in Fig. 7, nasal airflow continues almost throughout the following stop in descendent vers [desãd vER] 'going down to ' '7(ii) and in maisons qu'étaient [mEzO& ketE] 'house that were' (7(i)),even irrespective of a word boundary after the nasal Basset P., Amelot A., Vaissière J., Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100. vowel in the latter case. The maximal level of nasal airflow corresponds to the beginning of the oral closure. The burst is visible on the spectrogram and it still sounds like a stop (at least to the authors). As suggested by a rewiever of this paper, the upward excursion in the nasal graph in the /k/ case may be related to velic raising for stop closure, as the velum slides along the pharynx wall, pushed by tongue-dorsum/velum contact. This articulatory configuration acts like a piston creating a nasal air stream independent of air from the lungs. This may not be the only explanation, since the same phenomenon was observed for non velar There are some well-known cases of complete assimilation of unvoiced stops in French in words like "maintenant" [mE) tnA) ] 'now'. The pronunciation can vary from [mEt´nA) ] in hypoarticulated speech, to [mEnA) ] in spontaneous speech (figure 8). In the latter case, the voiceless stop is completely assimilated (or considered as absent). [START_REF] Dell | Deux nasalisations en français. Acte du séminaire « lexique et traitement automatique des langages[END_REF] postulated that no assimilation could occur in a sequence of [nasal] + [unvoiced stop] + [sonorant]. Our data (figure 7) provide counter evidence.

In N+V context, carryover seems to be the rule. (118/122, 97%).

Basset P., Amelot A., Vaissière J., Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100. Table 8 shows that there are no real differences between styles.

When compared with Cohn's results, our data differ only with regard to NAF before the nasal (see table 9). The present study could establish no clear tendency on either way. The results are in agreement for NAF spreading after the nasal.

Basset P., Amelot A., Vaissière J., Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100.

Table 9. Summary of the results obtained by [START_REF] Cohn | Phonetic and Phonological Rules of Nasalization[END_REF] in an aerodynamic study of words inserted in frame sentences in French compared to the results obtained in this study.

Contexts

Conclusion

It was hypothesised that less coarticulatory phenomena (anticipation and carry-over) would happen in read speech, as compared to spontaneous speech. Our data does not confirm such a trend.

The nasal airflow tends to appear before or at the onset of a nasal in both styles. A clear difference was observed depending on the identity of the surrounding phonemes. A majority of anticipatory phenomenon were found in the context C+Ṽ with a clear difference depending on the voicing of the consonant. When C is voiced there is a majority of occurrences where NAF precedes the nasal vowel. On the other end, a vast majority of items with NAF spreading after the offset of the nasal were observed in all contexts reported in this study. In conclusion, no clear differences could be established between speaking styles, neither for anticipatory nor carryover phenomena. More data are however needed to confirm the present results, to investigate more systematically inter and intra-speakers differences, the effect of rate of speech and styles, such as hyperarticulated and hypoarticulated speech. Perception tests are also needed to determine how nasalised vowels and consonants are perceived out of context. A reliable way remains to be found to compute the magnitude of nasal airflow in order to confirm [START_REF] Benguerel | Velar co articulation in French: an fiberscopic study[END_REF] finding that there is a difference in the opening of the velopharyngeal port between the syllable [na] and [nA) ].

Figure 1 .

 1 Figure 1. Waveform and aerodynamic data of the word "ton" [tO) ] 'your'. The arrows show that the stop burst is well aligned with the increasing of oral airflow. (Data collected with Eva workstation).

Figure 2 .

 2 Figure 2. Waveform and aerodynamic data for the word "ton" [tO) ] 'your'. The two arrows show that the stop burst is well aligned with the increase of oral airflow. (Data collected with PCquirer).

  Figures 4 and 5 illustrate anticipation and carryover nasalisation, respectively, with examples from spontaneous speech.

Figure 3 .

 3 Figure 3. Criteria for classifying nasal airflow as Anticipatory (a), synchronous, Delayed (d), and Carryover (c) based on the timing of the nasal airflow onset and offset relatively to the nasal phone boundary. (Dn = duration of the nasal phone; Dn-1 refers to the duration of the phone preceding the nasal and Dn+1refers to the duration of the phone following the nasal phone)

Figure 4 .

 4 Figure 4. Data illustrating the criteria used for characterizing anticipation. (i) synchronous onset of NAF, in the word "attention" [atA) sjO) ] 'be careful'. (ii) anticipation in "Provence"[pÂOvA) s]. (iii) delay in "d'embouts" [dA) bu] 'nozzles'. Data taken from spontaneous speech.

Figure 5 .

 5 Figure 5. Data illustrating the criteria used for characterizing carryover. (i) synchronous offset of NAF, in the word "bon" [bO) ) ] 'good'. (ii) carryover in "Provence"[pÂOvA) s]. Data taken from spontaneous speech.

  Basset P., Amelot A., Vaissière J., Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100.

Figure 7 .

 7 Figure 7. Example of (i)"maisons qu'étaient" [mEzO& ketE] 'houses that were'. (ii) "descendent vers" [desA) d vEÂ] 'going down to'.

Figure 8 .

 8 Figure 8. Example of the production of the same word "maintenant" [ mE& tnA& ] 'now' by the same speaker. The item on the left was extracted from a read corpus and the item on the right from spontaneous speech.

Table 1 :

 1 Calculations.

Table 2 : Results. Anticipation: Number of items, Mean a , ratio between the duration of anticipation and the total length of the preceding phone. Synchronous: number of items with NAF starting at the onset. Delay: number of items, Mean duration of delay, ratio of the duration of the delay on the total duration of the nasal

 2 Table5reveals no clear differences between read and spontaneous stylesBasset P., Amelot A., Vaissière J., Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100.Basset P., Amelot A., Vaissière J., Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100.

	.
				Anticipation			No anticipation
				NAF before onset		NAF	NAF after onset
	Contexts		Style	number of items	mean a	Dn a	-	1	Synchronous Number of items	number of items	Delay mean d	Dn d
			Spon	9	112ms	1.3		1	3	34 ms 0.3
		High vowels
	V+Ṽ	Read Low vowels Spon Read	8 -1	91ms --	0.9 --		4 2 1	1 1 1	33 ms 0.5 23 ms 0.4 18 ms 0.2
		Total		18					8	6
		Stops	Spon Read	14 9	71 ms 75 ms	1.1 1		1 2	1 5	47 ms 0.4 53 ms 0.3
	C voised +Ṽ	Fricatives	Spon Read	22 29	53 ms 70 ms	0.9 1		7 2	2 -	38 ms 0.3 --
		Laterals	Spon Read	2 2	76 ms 36 ms	2.7 1.1		1 -	-1	-64 ms 0.6 -
		Total		78					13	9
	C unvoiced +Ṽ	Stops Fricatives	Spon Read Spon Read	10 4 8 9	61 ms 16 ms 68 ms 77 ms	0.6 O.2 0.7 0.7	15 13 3 1	7 15 2 3	23 ms 0.2 45 ms 0.3 35 ms 0.2 46 ms 0.4
		Total		31					32	27
		High vowels Spon Read	24 29	77 ms 55 ms	1.2 0.8		33 26	4 6	32 ms 0.3 24 ms 0.4
	V+N	Low vowels Spon	9	54 ms	0.7		10	1	13 ms 0.1
			Read	13	59 ms	0.5		6	1	22 ms 0.4
		Total		75					75	12
				202					128	54

Table 3 :

 3 Total number of items in all contexts where NAF starts before, at the onset or after the onset of the nasal phone.

	Contexts Total number of items anticipation	synchronous	delay
			%	%	%
	V+Ṽ	32	56%	25%	19%
	C voised +Ṽ100	78%	13%	9%
	C unvoiced +Ṽ	90	34%	36%	30%
	V+N	162	46%	46%	8%
	Total	384	53%	33%	14%

Figure 6. Examples of anticipatory nasal airflow. (i) "C'est bon" [sE bO) ] 'is it all right'. (ii) "vacances" [vakA) s] 'holidays' (spontaneous speech)

Table 4 :

 4 Comparison between anticipation and no anticipation (both spontaneous and read speech styles merged).Basset P., Amelot A., Vaissière J.,Roubeau B. (2001), Nasal flow in French Spontaneous Speech, Journal of the International Phonetic Association, 31 (1) p. 87-100.

		Spontaneous		Read
		anticipation No anticipation anticipation No anticipation
	V+Ṽ	9	7	9	7
	C voised +Ṽ38	12	40	10
	C unvoiced +Ṽ	18	27	13	32
	V+N	33	48	42	39
	Total	98 (51%)	94(49%)	104(54%)	88(46%)

Table 6 :

 6 Results. No carryover: number of items. Carryover: number of items, Mean duration of carryover, ratio between the duration of carryover and the total length of the following phone.

				no carryover		carryover
							c
	Contexts		Style	nb of items	nb of items	Mean c	Dn	+	1
		High	Spon	5	2	149 ms	1.3
		vowels	Read	2	5	121 ms	0.8
	V+V	Low vowels	Spon Read	-2	4 2	110 ms 110 ms	1.3 0.1
		Total	9	13
		Stops	Spon Read	3 1	14 16	83 ms 88 ms		1 1
	V+C voiced	Fricatives Spon Read	2 -	12 14	62 ms 57 ms	1 0.7
		Laterals	Spon Read	--	4 4	130 ms 119 ms		2 2
		Total	6	64
	V+C unvoiced	Spon Read Fricatives Spon Stops Read	5 4 3 2	48 52 29 30	67 ms 56 ms 72 ms 69 ms	0.7 0.6 0.8 0.6
		Total	14	159
		High	Spon	2	50	121 ms	1.4
		vowels	Read	-	52	133ms	1.5
	N+V	Low	Spon	1	8	171 ms	1.7
		vowels	Read	1	8	121 ms	1.6
		Total	4	118
				33	354

Table 7 :

 7 Total number of items in all contexts where NAF ends at the offset of the nasal (no carryover) or spread on the following phone (carryover).

	Contexts	Total number of items No carryover	carryover
			%	%
	Ṽ +V	22	40%	60%
	Ṽ +C voised	70	8%	92%
	V+ C unvoiced	173	8%	92%
	N+ V	122	3%	97%
	Total	387	9%	91%

Table 8 :

 8 Comparison between styles.

		Spontaneous	Read
		No carryover carryover No carryover carryover
	V+V	5	6	4	7
	V+C	13	107	7	116
	N+V	3	58	1	60
	Total	21(11%)	171(89%)	12(6%)	183(94%)

Acknowledgements

We would like to thank J. J. Ohala for valuable comments. We wish also to thanks K. Kohler, A. Simpson and the reviewers for their helpful suggestions.


Xavier Dolques
email: xavier.dolques@engees.unistra.fr

Florence Le Ber
email: flo-rence.leber@engees.unistra.fr

Marianne Huchard
email: marianne.huchard@lirmm.fr

Clémentine Nebut
email: clementine.nebut@lirmm.fr

Relational Concept Analysis for Relational Data Exploration

Relational Concept Analysis (RCA) is an extension to the Formal Concept Analysis (FCA) which is an unsupervised classification method producing concept lattices. In addition RCA considers relations between objects from different contexts that allow for the creation of links between lattices. This feature makes it more intuitive to extract knowledge from relational data and gives richer results. However, data with many relations imply scalability problems and results that are difficult to exploit. We propose in this article a possible adaptation of RCA to explore relations in a supervised way in order to increase the performance and the pertinence of the results.

Introduction

Formal Concept Analysis [START_REF] Ganter | Formal Concept Analysis[END_REF]], written shortly FCA, is an automatic classification method of objects described by attribute through a binary relation. Such a classification results in a concept lattice (also called Galois lattice [START_REF] Barbut | Ordre et Classification: Algèbre et Combinatoire[END_REF]]) where each concept groups all the objects sharing the same attribute set. It is possible to navigate through a lattice in a simple and intuitive way, from the most specific concepts (concepts grouping many character-istics shared by only a few objects) to the less specific ones (concepts grouping many objects but sharing only a few characteristics).

FCA is used in several domains as a knowledge extraction method and the different publications on the topic, namely [START_REF] Carpineto | Concept Data Analysis: Theory and Applications[END_REF], [START_REF] Valtchev | Formal concept analysis for knowledge and data discovery : New challenges[END_REF] Valtchev, Missaoui, and Godin] have identified its forces and limitations. Some of those limitations can be worked around by using different approaches.

Relational Concept Analysis (RCA) [START_REF] Huchard | Relational concept discovery in structured datasets[END_REF]Huchard, Hacène, Roume, and Valtchev] is an extension of FCA taking into account relations between objects in addition to the characteristics of the objects. RCA consists in iteratively applying an FCA algorithm to deal with relational data: objects are described by attributes and by their relations towards other objects. Concepts discovered by a given iteration are propagated along the relations, leading to the discovery of new concepts at the next iteration.

RCA appears to be more intuitive to use on relational data such as databases or object-oriented modeling languages such as UML. In this article we propose to adapt RCA for the purpose of using it as a knowledge extraction method on water quality measures data for Alsacian watercourses.

This work is part of the ANR project FRESQUEAU1 which goal is to develop new study, comparison and exploitation approaches of all the available parameters on watercourses. It extends a previous study using FCA [Bertaux et al(2009a) Bertaux, Le Ber, Braud, and Trémolières] and statistical approaches [Bertaux et al(2009b) Bertaux, Le Ber, Li, and Trémolières].

Propagating along relations the concepts discovered from an iteration to another permits to discover interesting concepts, but it often leads to a combinatorial explosion of the number of concepts, and the interesting patterns are difficult to extract from the big concept set built. Several strategies can be used to counter this complexity, including the separation of the initial objects in several subsets after preliminary analysis or the introduction of requests [START_REF] Azmeh | Querying relational concept lattices[END_REF]. We are interested in this article in using RCA to explore interactively the data by letting the user choose before each iteration of FCA which contexts (object-attribute and object-object) he or she wants to use.

We are working on data that are not initially shaped as a binary relation but many works about data scaling will permit to get a binary relation [START_REF] Ganter | Formal Concept Analysis[END_REF]] or pattern structures [START_REF] Ganter | Pattern structures and their projections[END_REF]]. Those approaches have been previously applied on similar data in [Bertaux et al(2009a)Bertaux, Le Ber, Braud, and Trémolières] therefore in the following we only consider data as binary relations.

In this paper we present FCA, then the general principle of the RCA process as to highlight several variation points that would permit to improve the use of RCA in a data mining context. We then present an example of the kind of data from the FRESQUEAU project and the consequences of the variations on these data. We then conclude with a short discussion.

Formal Concept Analysis

FCA's purpose as defined by [START_REF] Ganter | Formal Concept Analysis[END_REF]] is to classify a set of objects described by attributes and presented as a formal context. A formal context K is a triplet (O, A, I) where O is an object set, A is an attribute set and I ⊆ O × A is the incident relation between O and A such that (o, a) ∈ I if and only if a is an attribute of o. Table 1 represent a formal context. The object set is here a set of identifiers for sampling stations on different watercourses. Each station is represented by a row. The attributes are description characteristics of the watercourses. The relation between a station and a characteristic of its watercourse is represented by a cross. Thus the station identified by BREI0001 is located in a small watercourse which water is fresh and live. The stations BRUN001 and BRUN002 are located on the same river but at different locations. Applying FCA on a context K = (O, A, I) leads to the generation of concepts. A concept is a couple (X,Y ) where X ⊆ O and Y ⊆ A such that X = {o ∈ O|∀a ∈ Y, (o, a) ∈ I} and Y = {a ∈ A|∀o ∈ X, (o, a) ∈ I}. X is called the extent of the concept Y its intent. The extent of a concept is the maximal set of objects sharing the intent attributes and the intent of a concept is the maximal set of attributes shared by all the extent objects.

For a given context FCA leads to the generation of all the concepts. Those concepts are forming a concept lattice also called Galois lattice. A concept c 1 is more general (resp. more specific) than a concept c 2 if the extent of c 1 contains (resp. is contained by) the extent of c 2 . In a dual way, the intent of a concept is contained by the intent of a concept more specific. Two given concepts have a unique superior bound and a unique inferior bound.

Lattices are usually represented by their Hasse diagram. The lattice of table 1 is represented by figure 1. Arrows are representing the generalization relation, i.e. the pointed concept is more general than the concept of the origin. Considering that the intent of a concept is included by the intent of every concept more specific and the extent of a concept is included by the extent of every concept more general, each object (resp. attribute) is displayed only once in the most specific (resp. most general) concept where it appears. For instance, concept 6 groups stations BREI0001 and DOLL001 that own attributes small watercourse and live and fresh water that can be found by the generalization relation towards concepts 2 and 7.

RCA Extension for exploratory analysis

Relational Concept Analysis (RCA) [START_REF] Huchard | Relational concept discovery in structured datasets[END_REF] is an extension of FCA considering, in addition to object characteristics, existing relations between objects.

The algorithm 1 present the main steps of RCA. The input parameter for RCA is a Relational Context Family RCF = (K, R) composed of n object-attribute contexts

K i = (O i , A i , I i ), i ∈ [1..n], and m object-object contexts R j = (O k , O l , I j ), j ∈ [1..m]
where O k and O l are object sets of K k et K l . It can be seen in table 2 an example of a Relational Context Family. It can be seen on the left hand side two objectattribute contexts taxons and stations and on the right hand side the objectobject context taxonPresence that links objects from context stations to objects of context taxons2 .

For R j ⊆ O k × O l , we call O k its domain and O l its range. The initialization step (lines 4-5) consists in building, for all i ∈ [1..n], the lattice L 0 [i] associated to the context K i . Figure 2 present the two lattices obtained after the initialization step on our example. It can be noticed that the relation taxonsPresence is not considered at this point of the processus and that the two lattices are independent.

At step p:

• EXTEND-REL add to K i the relations obtained by scaling the relations where K i is the domain. The scaling consists in the inclusion of the object-object rela-

1: proc MU L T I-FC A( Input: (K, R) a RCF, 2: Out: L table [1..n] of lattices) 3: p ← 0 ; halt ← false 4: for i from 1 to n do 5: L 0 [i] ← BU I L D-LA T T I C E(K 0 i ) 6:
while not halt do 7:

p + + 8: for i from 1 to n do 9: tions as relational attributes. They are obtained by using lattices concepts from step p -1 and a scaling operator (i.e. ∃, ∀). For example, if the scaling operator ∃ is chosen to scale a given relation R j , the columns of R j are replaced by attributes of the form ∃R j : C, where C is a concept in the lattice built from objects of the range of R j at step p -1. An object o from the domain of R j owns the relational attribute (∃R j : C) if R j (o) ∩ Extension(C) = / 0. Thus, we have extended the object-attribute context stations with the object-object context taxonsPresence scaled by the operator ∃. The extended context is presented by the table 3. The station FECH001 is linked to concepts 0, 1 and 3 by the relation ∃taxonsPresence as we can found the presence of Athericidae that can be found in concepts 0 and 1 and the presence of Boreobdella that can be found in concepts 0 and 3. If we had used the scaling operator ∀, the station FECH001 would be linked only to concept 0 by the relation ∀taxonsPresence as it is the only concept where we can find both taxons for this station.

K p i ← EX T E N D-RE L(K p-1 i , L p-1 ) 10: L p [i] ← UP D A T E-LA T T I C E S(K p i ,L p-1 [i]) 11: arrłt ← i=1,n IS O M O R P H I C(L p [i], L p-1 [i])
• UPDATE-LATTICE update the lattices of step p -1 to generate, for i ∈ [1..n], the lattice L p [i], associated to K i concatenated to every scaled object-object context which domain is K i .

The algorithm stops when a fix point is reached, i.e. when the lattice family obtained is isomorphic to the one from the previous step and the context extensions is unchanged. In our example, the lattices from figure 3 are the final lattices. The number of iterations is predictable when relations between contexts are not forming a circuit. But in some cases, for instance when an object-object context has same domain and range, the number of iterations is not predictable (only a maximal bound can be known) and can be really big depending on the data.

The relational lattices interpretation is different from the interpretation of classical concept lattices as they must be considered silmutaneously. The lattice stations of figure 3 must be considered with the lattice taxons to be correctly interpreted. We find in concept intents some attributes referring to other concepts. E.g. Concept 8 owns the relational attribute taxonsPresence : Concept 2 which means that all the objects of Concept 8 are linked by the relation taxonsPresence with at least (as the scaling operator used is ∃) an object of Concept 2 extent from lattice taxons.

The advantage of such a process is that the concepts obtained have in their intent relations to other concepts in addition to classical attributes. Those relations permit the extraction of patterns built from several interconnected contexts, as it has previously been done in [START_REF] Dolques | From transformation traces to transformation rules: Assisting model driven engineering approach with formal concept analysis[END_REF] Dolques, Huchard, and Nebut] and [START_REF] Dolques | Fixing generalization defects in UML use case diagrams[END_REF] Dolques, Huchard, Nebut, and Reitz], that could not be easily obtained from the classical FCA process.

However a major drawback of this kind of process is the potential difficulty to apprehend the result. In previous works in Model Driven Engineering, the data extracted from models of medium size can easily be apprehended by RCA. However, in a data mining context, the data size is more important. Computing time depends on the number of generated concepts, and it is exponential in regard to the minimum between the number of attributes or objects in the worst cas. Thus, if the relations between objects are numerous and with few similarities between objects, computation time can exponentially increase and the result may appear difficult to understand by a user because of the number of concepts to consider simultaneously. This is particularly true when only small patterns are needed when a lot of relations link the objects and these relations are forming a circuit. In such cases, we think it would be relevant to use and exploratory approach.

We list in the following the different possible variations on the algorithm to put in practice an exploratory approach. We enumerate the possible variation points in the algorithm that can affect the result by changing the contexts considered at each step. We propose for each variation point an alternate scenario from the previously described process that involve the user by asking him to choose. All those variations or only a subset of them can be applied depending on the needed granularity.

• initialization step, line 4 to 5 Build lattices for selected object-attribute contexts concatened to selected object-object contexts. • EXTEND-REL, line 9 Instead of using all the relations and scaling all the objectobject relation at each step, select a subset of the Relational Context Family and different scaling operators for each object-object context selected. Notice : the object-object relations need a lattice classifying the object of the range that must have been computed in a previous step, not necessarily p -1. At this step, objectattribute context can also be selected and the corresponding lattice can be built. In this section we illustrate the potential of an exploratory approach using the example of alsacian watercourses. Each watercourse is classified depending on general types (e.g. un small watercourse with fresh and live water). To evaluate the water quality of a watercourse, hydro-ecologists select several sections, called stations, on which they take samples and measures in addition to the sampling of the present plants and animals (called taxons). Those samples and measures respect several norms. After analysis or determination in a laboratory, the watercourse stations are described by different quantitative attributes : on one hand biological data (e.g. number of individuals for each taxon) and on the other hand physico-chemical data (e.g. pH, temperature, level of organic matter, level of dissolved oxygen, etc.). Taxons are characterized life traits that are qualitative data (e.g. lifetime of invertebrates). Figure 4 presents data with a schema. Each node represents an object set which is translated in an object-attribute context. Each labeled edge represents a relation between object sets which is translated in an object-object context. For each edge we consider the relation in both direction.

We would like to extract from these data some relations between the different kind of information that describe a station. It would be interesting for instance to extract rules with the format: modality M of life trait T in a watercourse of type E implies the presence of physico-chemical character C, which can be obtained by orienting thre relations of figure 4 such that the relation taxons' traits go from taxons to life traits and that all other relations have stations for domain as in figure 5.

Contextes objet-attribut

Station

ID BREI0001 ID BRUMB001 ID BRUN001 ID BRUN002 ID DOLL001 ID FECH001 BREI0001 × BRUMB001 × BRUN001 × BRUN002 × DOLL001 × FECH001 × traitDeVie ID inf-eg1 ID sup1 inf-eg1 × sup1 × taxons ID Asellus ID Athericidae ID Baetis ID Bithynia ID Boreobdella Asellus × Athericidae × Baetis × Bithynia × Boreobdella × descriptionCours ID petit
× caracPC ID DCO ID NH4 ID MES ID T DCO × NH4 × MES × T × Contextes objet-objet caracPC1 DCO NH4 MES T BREI0001 × × × × BRUMB001 × × BRUN001 × × BRUN002 × × × DOLL001 × × FECH001 × × caracPC2 DCO NH4 MES T BREI0001 BRUMB001 × × BRUN001 × × BRUN002 × DOLL001 × × FECH001 × × taxon-trait inf-eg1 sup1 Asellus × Athericidae × Baetis × Bithynia × × Boreobdella × station-descCours petit cours d'eau grand cours d'eau eau calme et frache eau calme et vive cours phratique BREI0001 × × BRUMB001 × × BRUN001 × BRUN002 × DOLL001 × × FECH001 × × station-taxon-20-99 Asellus Athericidae Baetis Bithynia Boreobdella BREI0001 BRUMB001 × BRUN001 BRUN002 × DOLL001 × FECH001 station-taxon-1-19 Asellus Athericidae Baetis Bithynia Boreobdella BREI0001 × × × BRUMB001 BRUN001 × × × × BRUN002 × × × DOLL001 × × FECH001 × × × × station-taxon-100+ Asellus Athericidae Baetis Bithynia Boreobdella BREI0001 × × BRUMB001 BRUN001 BRUN002 DOLL001 FECH001
Table 4 Famille relationnelle de contextes obtenue partir de nos donnes d'exemples en considrant les directions de relation de la figure 5.

The Relational Context Family for this particular configuration is presented by table 4. We add in the object-attribute contexts an identifier as attribute for each object in order to create a concept for each object after the initialization step. Figure 6 presents needed excerpts of the lattice family obtained to extract the rule Presence at medium level of a taxon which lifetime is more than one year implies a high Chemical oxygen demand. 3 Indeed, we can notice that Concept 3 is more specific that Concept 41. An attribute introduced by Concept 3 implies an attribute introduced by Concept 41 as all the attributes of Concept 41 are inherited by Concept 3. Concept 3 groups the stations that shelter a medium quantity of individuals from taxons grouped by Concept 63. Concept 63 groups taxons which individuals have a lifetime superior to one year. Concept 41 groups the stations which physico-chemical characters (of level 2, i.e. a high level) is Chemical oxygen demand. From those observations we can induce the previous implication. From the complete lattice we can obtain the whole set of implication rules between life traits and physico-chemical characters by considering all the cases where physico-chemical characters are introduced by a concept and life traits are introduced by a more specific concept.

But rules in the following format can also be relevant: modality M of life trait T can appear when the physico-chemical character C is present. To generate these kind of rule, we need to change from the previous configuration the direction of the relations between life traits and taxons and between taxons and stations. There exists more configurations which results can be relevant and varying the scaling operators increases the expression of the rules obtained with RCA. If we consider the schema from figure 4 as a graph, the exploration consists in analyzing the different edges until we obtain a path between the physico-chemical characters and the life traits. By comparing the obtained results during an exploration with the results obtained using the classic RCA on the same relation, we notice that the final lattices are smaller and easier to read in the first case.

Combining the whole set of possible configurations is not to consider as by multiplying the relations we take the risk to raise a combinatorial explosion of the number of concepts, thus increasing the computation time and the complexity of the obtained concepts. For our example, considering all the relations in both directions scaled only by ∃ lead to the creation of 120 concepts against 66 and 63 concepts for the two configuration previously presented. So we plan an approach where the user can explore different configuration by making different choices at each step of the process as presented in section 1.

Conclusion and discussion

In this article, we have presented an exploratory approach to assist the use of RCA in a way more appropriated for knowledge retrieval process. We have several motives to modify the original RCA process: to obtain relevant results faster by computing less lattices (preferably only the lattices that are of interest), to decrease the complexity of relational data mining, or to let the expert guide the discovery process based on its intuition and the learning patterns that appear along the process.

Several questions arise on this concept extraction approach from relational data. The initialization step strongly impact the discovered structures. It can speed up the process, if the object-object relations contain the information needed by the expert, or on the contrary it can hide relevant information to the expert. Nevertheless, the most important problem comes from the fact that modifications at each step implies that the concept generation is not monotonous anymore and it becomes possible to build examples where the process diverges by iterating on several recurring configurations.

In the original RCA process, when the fix point is reached the lattices on the two last step are isomorphic. So when a concept refers another one via a relational attribute the referenced concept can be found in a lattice from the same step. But with the exploratory process that we are proposing, when a concept refers to another, the referenced concept is in a lattice of the previous step and this concept can also refers to a concept from a lattice in another previous step. So we need to find solutions to present the information to the expert in a way simple enough to interpret. However, we think that such an exploratory approach is more applicable than a systematic approach that iterate until a fix point and give results too difficult for an expert to interpret.

Fig. 1

 1 Fig.1Hasse diagram of the concept lattice generated from the context described by table 1.

Fig. 4

 4 Fig.4Schema of the data used by our example.

Fig. 5

 5 Fig. 5 Schema of data used by RCA.

Table 1

 1 Example of formal context. Objects are presented as rows and attributes as columns.

		small	large	fresh and calm fresh and live phreatic
		watercourse watercourse	water	water	watercourse
	BREI0001	x			x
	BRUMB001	x		x
	BRUN001				x
	BRUN002				x
	DOLL001	x			x
	FECH001		x		x

Table 2

 2 Algorithme 1: Processus of Relational Concept Analysis. Relation Context Family example. Objects are presented as rows and attributes as columns.

		object-attribute contexts	object-object contexts
	taxons ≤1 year > 1 year Athericidae x Bithynia x x Boreobdella x	taxonPresence BREI0001 BRUN001 FECH001	Atheri-Bithy-Boreob--cidae -nia -della x x x x x
	stations	small watercourse watercourse watercourse fresh and live phreatic
	BREI0001	x	x
	BRUN001			x
	FECH001		x
				Concept_4
		Concept_0	Concept_6	Concept_8
				fresh and live water	phreatic watercourse
				FECH001	BRUN001
		Concept_1	Concept_3	Concept_5
		≥ 1 year	< 1 year	small watercourse
		Athericidae	Boreobdella	BREI0001
		Concept_2	Concept_7
		Bithynia
		taxons	stations
	Fig. 2 Lattice Family generated from the Relational Context Family of table 2 after initialization
	step 0.

Table 3

 3 Scaling of the relation taxonPresence and extension of context stations at step 1. Lattice Family generated from the Relational Context Family of table 2 after step 1.

				Concept_4
			∃ taxonPresence : Concept_0
			∃ taxonPresence : Concept_3
	Concept_0	Concept_6	Concept_10
			fresh and live water	∃ taxonPresence : Concept_1
	Concept_1	Concept_3	Concept_5	Concept_9	Concept_8
	≥ 1 year	< 1 year	small watercourse	phreatic watercourse ∃ taxonPresence : Concept_2
	Athericidae	Boreobdella	BREI0001	FECH001	BRUN001
	Concept_2		Concept_7
	Bithynia
	taxons		stations
	Fig. 3

  • UPDATE-LATTICES, line 10 Update only the lattices for the selected relations.• stop, line 11 If a fix point is not reached, let the stop decision to the expert.

	4 Exploration example
						présence de
						caractères
						physico-chimiques
						niveau 1
		traits des		presence de taxons 1-19 presence de		présence de caractères physico-chimiques	caractères physico-chimiques
	traits de vie	taxons	taxons	taxons 20-99 presence de	stations	niveau 2 type de
				taxons 100-+		cours d'eau	types de
							cours d'eau

  Excerpts of lattices obtained from the relational context family of table 4.

	Extrait du treillis du contexte
	stations
	Extrait du treillis du contexte
	types de cours d'eau
		Extrait du treillis du contexte
	Extrait du treillis du contexte	caractères physico-chimiques
	taxons
	Extrait du treillis du contexte
	traits de vie
	Fig. 6

http://engees-fresqueau.unistra.fr/

The term taxon covers diverse terms used for the denomination of living beings such as species, genus or families.

considering the size of the example, this rule is to be considered for illustrative purpose only.


Miguel Angel Veganzones

Sylvain Douté

Jérémy E Cohen

Rodrigo Cabral Farias

Jocelyn Chanussot

Pierre Comon

R Cabral Farias

Nonnegative CP decomposition of multiangle hyperspectral data: a case study on CRISM observations of martian icy surface

Keywords: CRISM mission, multi-angle hyperspectral, MARS-ReCO, Canonical-Polyadic, nonnegative tensor decomposition

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

INTRODUCTION

The Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) sensor aboard the Mars Reconnaissance Orbiter takes hyperspectral acquisitions of Martian surface from the top of the atmosphere (TOA) using 11 viewing angles in 437 spectral bands, on visible and infrared wavelengths. The CRISM experiment is an invaluable source of information related to the scattering of solar light by surface materials and atmospheric aerosols. However, the extraction of surface bidirectional reflectance (BRF) is especially challenging, due to the strong effects of aerosols in Martian atmosphere. The Multiangle Approach for Retrieval of Surface Reflectance from CRISM Observations (MARS-ReCO) [START_REF] Ceamanos | Surface reflectance of Mars observed by CRISM/MRO: 1. multi-angle approach for retrieval of surface reflectance from CRISM observations (MARS-ReCO)[END_REF][START_REF] Douté | Multi-angle approach for coherent retrieval of surface reflectance and aerosol optical depth from CRISM observations[END_REF] defined an innovative TOA radiance model and inversion scheme aimed at correcting for aerosols effects taking advantage of the near-simultaneous multi-angle CRISM observations. MARS-ReCO retrieves the aerosol optical depth (AOD) and the BRF of surface materials as a function of wavelength.

The validation of MARS-ReCO has been solely performed on the base of controlled synthetic data [START_REF] Ceamanos | Surface reflectance of Mars observed by CRISM/MRO: 1. multi-angle approach for retrieval of surface reflectance from CRISM observations (MARS-ReCO)[END_REF] and by comparison with in situ reflectance measurements by the Martian Exploration Rovers (MERs) [START_REF] Fernando | Surface reflectance of Mars observed by CRISM/MRO: 2. estimation of surface photometric properties in gusev crater and meridiani planum[END_REF]. This validation was restricted to moderately anisotropic mineral surfaces. Here, we aim to provide supporting evidence by performing a spectral unmixing (SU) analysis of MARS-ReCO BRF estimates from CRISM multiangle observations of highly reflective and anisotropic icy surfaces at high latitudes with grazing illumination. This is part of a collaboration effort to validate MARS-ReCO through the cross-validation of different techniques in the absence of ground truth [START_REF] Douté | Testing the multi-angle approach for retrieval of surface reflectance from crism observations on martian icy surfaces[END_REF]. Spectral unmixing (SU) is one of the most studied topics in hyperspectral image analysis. The goal of SU techniques is to estimate the spectral signatures of the materials present in the image (endmembers) and their spatial distributions (fractional abundances). Usually, a linear mixing model (LMM) is considered:

X = ΦS + E, (1)
where S ∈ R L×R + denotes the matrix of endmembers, each of the R columns representing the spectral signature of a given macroscopic material, Φ ∈ R N ×R + denotes the matrix of fractional abundances corresponding to each of the R endmembers; and E ∈ R N ×L denotes additive noise.

In [START_REF] Veganzones | Multilinear spectral unmixing of hyperspectral multiangle images[END_REF][START_REF] Veganzones | Nonnegative tensor CP decomposition of hyperspectral data[END_REF] we proposed to exploit the hyperspectral multiangle information by representing the hyperspectral multiangle acquisitions as a nonnegative tensor, X ∈ R N ×L×T + , with N pixels, L wavelengths and T angles of acquisition. Tensor analysis [START_REF] Comon | Tensors : A brief introduction[END_REF] (a.k.a. multiway or multiarray analysis) techniques, and specifically the Canonical Polyadic (CP) model [START_REF] Hitchcock | The expression of a tensor or a polyadic as a sum of products[END_REF], sometimes coined Candecomp/Parafac [START_REF] Kiers | Towards a standardized notation and terminology in multiway analysis[END_REF], can be used to decompose the multiangle tensor, X , into a multilinear combination of R nonnegative factor matrices A, B and C, of size N × R, L × R and T × R, respectively:

X ijk = R r=1 AirBjrC kr , (2)
where the A, B and C factors are associated with estimations of the fractional abundances, endmembers and angular anisotropies of the materials in the scene. Thus, the nonnegative CP decomposition can be understood as a multilinear spectral unmixing, exploiting the additional angular information.

The rank of tensor X is defined as the minimal number R of terms necessary for the equality above to hold exactly. A shorter expression can be employed to denote the decomposition in (2) using a diagonal tensor of ones, I, of size R×R×R, i.e., (A ⊗ B ⊗ C) I, where ⊗ is a tensor product. In practice, the data tensor is subject to modeling errors or measurement noise, and it is convenient to find its best rank-R approximation by minimizing an objective function of the form

Υ(A, B, C) = X -(A ⊗ B ⊗ C) I 2 , (3)
for some well chosen norm, instead of solving the exact CP decomposition [START_REF] Douté | Multi-angle approach for coherent retrieval of surface reflectance and aerosol optical depth from CRISM observations[END_REF]. It is known that tensors of order 3 or larger do not always admit a best rank-R approximate, when R > 1, especially in R or C. But fortunately, it has been shown in [START_REF] Lim | Nonnegative approximations of nonnegative tensors[END_REF] that such obstacle does no exist for nonnegative tensors, and that the problem is well-posed in R + : best lower rank nonnegative approximates always exist and are generically unique under mild conditions on R [START_REF] Qi | Uniqueness of non-negative tensor approximations[END_REF].

Here, we study the application of the nonnegative CP decomposition method to unmix a multiangle BRF data cube extracted by MARS-ReCO from a challenging area of interest in the outskirts of the south permanent polar cap in Mars.

DATA AND METHODS

Data

The data set is a multi-angle CRISM acquisition in the outskirts of the south permanent polar cap in Mars. The analysis here is restricted to the visible channel (VIS) of the instrument (0.4-1.0µm). The scene features regions of submiating icy CO2, revealing H2O frost and defrosted terrains. The topography of the region is moderate with maximum slopes of 10%. No ground truth is available for this region. Table 1 gives information about the data set, and Fig. 3 depicts a high-resolution gray scale image of the scene taken at nadir. In order to facilitate the simultaneous processing of the multiangle TOA reflectance factor of CRISM observation, the spectra corresponding to the eleven view angles are rearranged in a common geographical grid of super-pixels. Each spatial bin is approximately 200 meters in size [START_REF] Douté | Multi-angle approach for coherent retrieval of surface reflectance and aerosol optical depth from CRISM observations[END_REF]. Super pixels are stacked in one way, so that after preprocessing the resulting data is a third order tensor with super pixels, wavelengths and angles as diversities.

Nonnegative CP decomposition

Given a nonnegative tensor, X ∈ R N ×L×T + , its approximation by a nonnegative CP decomposition is formulated as:

argmin X -(A ⊗ B ⊗ C) I 2 F w.r.t. A, B, C s.t. A 0, B 0, C 0, (4)
where • F denotes the Frobenius norm, and inequality applies entry-wise. This problem is highly non-convex. The most commonly used method for nonnegative CP decomposition is alternating nonnegative least squares (ANLS) [START_REF] Cichocki | Nonnegative matrix and tensor factorizations[END_REF], where the cost function is minimized in an alternating way for each factor (A, B or C) while the others are fixed.

When the data tensor is large as in hyperspectral multiangle data, where the number of pixels is usually larger than 10 6 , the workhorse techniques described briefly above can fail to handle all the data within the memory of the computer, or can converge very slowly. An approach to handle large tensor decomposition is through the use of compression. The general idea is that the original tensor X can be equivalently represented by a compressed version of it, X c, with reduced dimensions Nc × Lc × Tc. The compressed tensor is then decomposed by minimizing

Υ = X c -(Ac ⊗ Bc ⊗ Cc) I 2 F w.r.t. Ac, Bc, Cc, (5)
where Ac, Bc, Cc are compressed versions of the original factor matrices. Note that, after the compressed factors are obtained, a decompression operation must carried out to recover nonnegative factors in the original dimensions.

Authors in [START_REF] Cohen | Fast decomposition of large nonnegative tensors[END_REF] provide two algorithms, a compressed conjugate gradient (CCG) and a Projected and Compressed Alternating Least Squares (ProCo-ALS) to solve (5) when the compression is done using an approximate High Order Singular Value Decomposition (HOSVD):

X ijk ≈ Nc,Lc,Tc lmn U il V jmW kn [X c] ijk , (6)
or more compactly X ≈ (U ⊗ V ⊗ W ) X c, where U , V and W are matrices with orthogonal unit-norm columns. These algorithms solve the optimization problem defined in (5) enforcing the nonnegativity of the factors in the unconstrained domain:

A ≈ U Ac 0, B ≈ V Bc 0, C ≈ W Cc 0. (7)
Here, we make use of a modified version of ProCo-ALS to include the sum-to-one physical constraint on abundances, usually employed in spectral unmixing:

A1R = 1N . (8)
The modified ProCo-ALS algorithm works by projecting A onto the unitary simplex instead of projecting it onto the nonnegative orthant using the solution to the waterfilling problem.

EXPERIMENTAL METHODOLOGY AND RESULTS

We run two sets of 20 Monte Carlo runs of the modified ProCo-ALS CP decomposition over the Mars reflectance tensor representation, using two tensor rank values, R = {3, 4}. Those superpixels of the reflectance tensor for which there were missing values were removed, so the number of tensor superpixels was N = 1633. The dimensions of the compressed tensor were set to Nc = 550, Lc = 50 and Tc = 11, so only the spatial way was slightly compressed in approximately a 1/3 ratio. For each ProCo-ALS run, we randomly initialized the factors A, B, C to nonnegative matrices, ensuring that the spectral factor, B, is normalized and, that the spatial factor, A, holds the sum-to-one constraint. Among the 20 runs, we kept the one maximizing the volume of the simplex formed by the columns of the spectral factor matrix, V (B) [START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[END_REF]:

V (B) = 1 R -1! det 1 . . . 1 b1 . . . bR . (9)
Figs. The VIS range offers no distinct features to distinguish the materials present on the image, because of the relative flat behaviour of the spectral factors. On the contrary, the angular factors give information that can be physically interpreted. The first two compo- nents present high anisotropic behaviours. The first one displays a strong backward scattering lobe (high values at low phase angles, T = 1 . . . 5), associated with the presence of extended shadows (sun low on the horizon) at the sub-pixel scale. The second one has a marked forward-scattering scattering lobe (high values at high phase angles, T = 8 . . . 11), linked with volume scattering within the surface materials. Finally, the third component shows much less anisotropic behaviour and, upon mixing, likely controls the overall level of the BRF as a function of geometry, i.e. the intrinsic albedo. Fig. 5 shows the results of a k-means clustering over the spatial factors obtained by ProCo-ALS for a number of clusters k = 5. The fourth component estimated by ProCo-ALS with R = 4 has been discarded in the analysis. Results obtained by setting k = 5 present a high degree of internal coherence..According to Fig. 5, the dark unit of the nadir image (top and middle-left strips) could be decomposed into two clusters (red and green), while the brightest unit of the nadir image consists in three distinguishable clusters (blue, cyan and yellow). The latter three clusters share the high presence of the third component where dominates highly reflective CO2 ice but with a moderately anisotropic behaviour: smaller relative proportions of the other two materials (with high anisotropies). The difference among these segmented areas appears to be related to the topography of the region and its relation to the shadowing effects. The dark unit, composed of two distinct areas, presents a higher relative abundance of the second component (backward scattering) but a varying abundance of the third component (moderate for the fourth cluster (red) but high for the fifth (green). Also, the former cluster presents a higher abundance of first component relative to the lat- ter cluster. The green cluster corresponds to a mixture of dust and minute amount of water ice (mostly defrosted terrains), while the red cluster corresponds to compositional transition zones with a high concentration of H2O frost. These results are in accordance with the ones obtained from a temporal series of the studied area in [START_REF] Douté | Testing the multi-angle approach for retrieval of surface reflectance from crism observations on martian icy surfaces[END_REF]. As expected, the results of the unmixing algorithm for the spectral way are hard to interpret since in the visible range, materials expected to be found on the surface of Mars present no discriminative features. Moreover in our analysis, the recovered spectra are almost collinear, suggesting the spectral diversity may not be necessary. Further analysis will be carried out to try other low rank unmixing methods like Nonnegative Matrix Factorization. Coupling present results with an analysis of infrared data for which materials should have distinguishable spectra should also prove interesting.

CONCLUSIONS

We have studied the unmixing of a multiangle BRF data cube extracted by MARS-ReCO from a challenging area of interest in the outskirts of the south permanent polar cap in Mars, by means of nonnegative CP decomposition. This work is part of a collaborative effort to validate MARS-ReCO through the cross-validation of different techniques in absence of ground truth. Obtained results are in accordance with other studies involved in the MARS-ReCO validation project, supporting the estimations of MARS-ReCO and the use of nonnegative CP decompositions to unmix multiangle hyperspectral data.

Fig. 3 .

 3 Fig. 3. High-resolution gray scale image of the study case scene taken by CRISM sensor at nadir angle.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Best run ontained by the modified ProCo-ALS with R = 3: (top) Spectral factor, B, (mid) Angle factor, Ĉ, and (bottom) Spatial factor, Â.

  Figs.1 and 2show the factors estimated by the best run according to[START_REF] Kiers | Towards a standardized notation and terminology in multiway analysis[END_REF] for R = 3 and R = 4, respectively. The only significant differences between the two solutions are: (i) the angle factors of the first and second components have a low tail appearing at last angle indices of the first component for R = 3, and at the first angle

Fig. 4 .

 4 Fig. 4. Results of a k-mean clustering (k = 5) over the spatial factor estimated by ProCo-ALS with (a) R = 3 and (b) R = 4. Figures depict the clustering maps in transparent colors over a high resolution grayscale image of the scene.

Fig. 5 .

 5 Fig. 5. Results of a k-mean (k = 5) over the spatial factor estimated by ProCo-ALS with (a) R = 3 and (b) R = 4. Bar plots show the centroids of the clusters in terms of the three components estimated by ProCo-ALS (the fourth component for R = 4 has been discarded in the analysis).

Table 1 .

 1 Information of the CRISM observation used in the study case.

	CRISM ID Ls 1	Incidence Phase range AOD 2
	FRT144E9 320 • 72 •	35-128 •	0.62

This work is supported by the European Research Council under the European Communitys Seventh Framework Programme FP7/2007-2013 Grant Agreement no. 320594, "DECODA". Also, by a contract with PNTS, ANR (AAP Blanc 2012 Inter II -SIMI 5) and a grant from Label OSUG@2020 (Investissements davenir ANR10 LABX56).


Miguel Angel Veganzones

Jérémy E Cohen

Rodrigo Cabral Farias

Konstantin Usevich

Lucas Drumetz

Jocelyn Chanussot

Pierre Comon

Rodrigo Cabral

R Cabral Farias

Canonical polyadic decomposition of hyperspectral patch tensors

Keywords: Spectral unmixing, extended linear mixing model, Canonical Polyadic, nonnegative tensor decomposition, patch tensor

INTRODUCTION

Imaging spectroscopy [START_REF] Goetz | Imaging spectrometry for earth remote sensing[END_REF] (a.k.a. hyperspectral imaging) is concerned with the measurement, analysis, and interpretation of spectra acquired from a given scene or object [START_REF] Plaza | Recent advances in techniques for hyperspectral image processing[END_REF]. Given a nonnegative hyperspectral image, X ∈ R N ×L + , where N denotes the number of pixels and L the number of spectral bands, spectral unmixing (SU) aims at estimating the spectral signatures of the materials present in the image and their spatial distributions, known respectively as endmembers and fractional abundances. This analysis ultimately permits a sub-pixel resolution. Usually, a linear mixing model (LMM) is considered:

X = ΦS + E, (1)
where S ∈ R L×R + denotes the matrix of endmembers, each of the R columns representing the spectral signature of a given macroscopic material, Φ ∈ R N ×R + denotes the matrix of fractional abundances corresponding to each of the R endmembers; and E ∈ R N ×L denotes additive noise. The endmembers are usually extracted from the data using one of the many endmember extraction algorithms available in the literature [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regressionbased approaches[END_REF]. Then, the abundances are estimated by solving an inverse problem. For instance, the full constrained least squares unmixing (FCLSU) is defined as:

argmin X -ΦS 2 F w.r.t. Φ s.t. Φ 0, Φ1 R = 1 N , (2)
with 1 K is an unitary vector of size K.

Sometimes, the LMM is insufficient, especially when the macroscopic materials features spectral variabilities due to illumination inhomogeneities, topography and other effects. Recently, the extended LMM (ELMM) [START_REF] Veganzones | A new extended linear mixing model to address spectral variability[END_REF] was proposed to address the spectral variability issue:

X = Φ f k (S) + E, (3)
The function f (•) is a mapping f : R L×R + → R L×R + that accounts for the spectral variability of the endmembers. The k index of the spectral variability function accounts for different flavours of the ELMM. For instance, the spectral variability function could vary across the pixels.

We propose to use tensor analysis [START_REF] Comon | Tensors : A brief introduction[END_REF] (a.k.a. multiway or multiarray analysis) to decompose the hyperspectral image in a few rank-one terms, by means of the Canonical Polyadic (CP) model [START_REF] Hitchcock | The expression of a tensor or a polyadic as a sum of products[END_REF], sometimes coined Candecomp/Parafac [START_REF] Kiers | Towards a standardized notation and terminology in multiway analysis[END_REF]. Recently, we have given evidence that the CP decomposition of hyperspectral time series and multiangle acquisitions could be understood as the multilinear spectral unmixing of the data [START_REF] Veganzones | Nonnegative tensor cp decomposition of hyperspectral data[END_REF]. However, it is not possible to apply this analysis directly to a single hyperspectral image. Some authors have proposed to consider the hyperspectral data cube as a tensor representation of the data, where the two spatial dimensions, i.e., rows and columns of the image, define two different ways of the tensor [START_REF] Zhang | Tensor discriminative locality alignment for hyperspectral image spectralspatial feature extraction[END_REF][START_REF] Bourennane | Improvement of classification for hyperspectral images based on tensor modeling[END_REF]. However, this representation is of little use for the CP decomposition since the resulting tensor is not low rank. Thus, we propose to build what we have termed a patch-tensor representation of a hyperspectral image by stacking small patches of the image. The resulting tensor, X ∈ R N ×L×B + , has a third way containing the B neighbouring pixels of each of the N pixels of the image. In Sec. 2, we provide a detailed explanation of the proposed patch-tensor representation.

The nonnegative CP decomposition allows to decompose the patch-tensor, X , into a multilinear combination of R non-negative factor matrices A, B and C, of size N × R, L × R and B × R, respectively, and :

X ijk = R r=1 A ir B jr C kr , (4)
where the scaling of factors has been pulled in C. The rank of tensor X is defined as the minimal number R of terms necessary for the equality above to hold exactly [START_REF] Comon | Tensors : A brief introduction[END_REF]. A shorter expression can be employed to denote decomposition (4) using a diagonal tensor of ones, denoted I, of size R × R × R, i.e. X = (A ⊗ B ⊗ C) I, where ⊗ is the tensor product.

In practice, the data tensor is subject to modeling errors or measurement noise, and it is convenient to find its best rank-R approximation by minimizing the following objective function

Υ(A, B, C) = X -(A ⊗ B ⊗ C) I , (5)
for some well chosen norm, instead of attempting to compute the exact CP decomposition of X with, unavoidably, a large rank. It is now known that tensors of order 3 or larger do not always admit a best rank-R approximate when R > 1, especially in R or C. But fortunately, it has been shown in [START_REF] Lim | Nonnegative approximations of nonnegative tensors[END_REF] that this obstacle does no longer hold for nonnegative tensors, and that the problem is well-posed in R + : best lower rank nonnegative approximates always exist and are generically unique under mild conditions on R [START_REF] Qi | Uniqueness of nonnegative tensor approximations[END_REF]. We show in Sec. 3, that the nonnegative CP decomposition of a hyperspectral patch-tensor is equivalent to solving a blind regularized version of the ELMM model [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regressionbased approaches[END_REF], where rank-one factor estimates correspond to spatial abundances, endmembers and local spectral variabilities. The CP decomposition is a blind technique in the sense that no a priori information is needed, i.e., the spectral signatures of the materials in the image are assumed unknown. We further discuss on some preliminary experimental results using synthetic data in Sec. 4, and we finally give some conclusions in Sec. 5.

PATCH-TENSOR REPRESENTATION

Let us denote by x i ∈ R L + , the i-th pixel, 1 ≤ i ≤ N , of a hyperspectral image. For each pixel x i , a patch is defined by the matrix

P i = [x i1 , x i2 , . . . , x i B ], of B pix- els x i k in the neighbourhood of x i , including x i = x i1 .
A neighbourhood of x i can be the set of pixels adjacent to x i using a sliding window (see Fig. 1). Each matrix

P i ∈ R L×B + is stacked to define a slice of the patch-tensor, X ∈ R N ×L×B + . The above formulation of the patch- tensor is equivalent to stacking B images X (k) ∈ R N ×L + , 1 ≤ k ≤ B,
obtained by shifting the original hyperspectral image along the neighbourhood, i.e. a shifted image is defined by X (k) = x (1+d k ) , . . . , x (N +d k ) , where d k denotes a spatial displacement from the center of the patch to the position of the k-th element in the patch. This yields zero-padding along the border of the hyperspectral image.

x i1 x i2 x i3 x i4 x i5 x i6 x i7 x i8 x i9 x i10 x i11 x i12 x i13 x i14 x i15 x i16 x i17 x i18 x i19 x i20 x i21 x i22 x i23 x i24 x i25 L x i1 x i2 x i3 x ik x i24 x i25 1 B 1 L P i Fig. 1.
Example of a 5 × 5 patch using a sliding window centred in x i , and the corresponsing matrix patch, P i .

NONNEGATIVE CP DECOMPOSITION OF THE PATCH-TENSOR 3.1. Compressed nonnegative CP decomposition

Given the patch-tensor, X ∈ R N ×L×B + , its approximation by a nonnegative CP decomposition ( 4) is formulated as:

argmin X -(A ⊗ B ⊗ C) I 2 F w.r.t. A, B, C s.t. A 0, B 0, C 0, (6)
where • F denotes the Frobenius norm. This problem is highly non-convex, yet many algorithms provide rather precise but costly computation.

An approach to handle large tensor decomposition is through the use of compression. The general idea is that the original patch-tensor X can be equivalently represented by a compressed version of it, X c , with reduced dimensions N c × D c × T c . The compressed tensor is then decomposed by solving

argmin Υ = X c -(A c ⊗ B c ⊗ C c ) I 2 F w.r.t. A c , B c , C c , (7)
where A c , B c , C c are compressed versions of the original factor matrices. Note that, after the compressed factors are obtained, a decompression operation is carried out to recover the factors in the original dimensions.

Authors in [START_REF] Cohen | Fast decomposition of large nonnegative tensors[END_REF] provide two algorithms, a compressed conjugate gradient (CCG) and a Projected and Compressed Alternate Least Squares (ProCo-ALS) to solve [START_REF] Kiers | Towards a standardized notation and terminology in multiway analysis[END_REF] when the compression is done using an approximate High Order Singular Value Decomposition (HOSVD):

X ijk ≈ Nc,Dc,Tc lmn U il V jm W kn [X c ] ijk , (8)
or using the same notation as in the CP model

X ≈ (U ⊗ V ⊗ W ) X c , (9)
where U , V and W are matrices with orthogonal unit-norm columns. These algorithms solve the optimization problem defined in [START_REF] Kiers | Towards a standardized notation and terminology in multiway analysis[END_REF] enforcing the nonnegativity of the factors in the unconstrained domain:

A ≈ U A c 0, B ≈ V B c 0, C ≈ W C c 0. (10
) Here, we make use of a modified version of ProCo-ALS to include the abundances sum-to-one physical constraint usually employed in spectral unmixing:

A1 R = 1 N , (11)
as in [START_REF] Plaza | Recent advances in techniques for hyperspectral image processing[END_REF]. The modified ProCo-ALS algorithm works by projecting onto the unitary simplex instead of projecting onto the nonnegative orthant.

Regularized extended linear mixing model (rELMM) equivalence

Here, we show that the CP decomposition ( 6) of the patchtensor, X , is equivalent to a regularized version of an ELMM formulation [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regressionbased approaches[END_REF]. Given the tensor slice, X :,:,k ≡ X (k) , where the colon symbol denotes all the elements of the corresponding way, its CP decomposition is given by:

X :,:,k ≈ R r=1 A :,r B :,r C k,r = AΛ k B , (12)
Λ k being a R × R diagonal matrix with entries diag (Λ k ) = C k,: . If we consider that factor matrices A and B are respectively estimates of the abundance matrix Φ and the set of endmembers S, then the CP model in ( 12) is equivalent to:

X (k) ≈ ΦΛ k Ŝ . (13)
Thus, the nonnegative CP decomposition of the patchtensor X can be reformulated as a blind joint solution to the following B ELMM models:

X (k) ≈ Φ f k (S) , k = 1, . . . , B, (14)
where the spectral variability function varies for each of the B slices:

f k (S) = SΛ k , k = 1, . . . , B. (15)
That is, the CP decomposition provides a blind solution to the following optimization problem:

argmin X (1) -Φ f 1 (S) 2 F + . . . + X (B) -Φ f B (S) 2 F w.r.t. Φ, S, Λ 1 , . . . , Λ B s.t. Φ 0, Φ1 R = 1 N , S 0, Λ 1 0, . . . , Λ B 0.
(16) The model in (16) corresponds to a regularized ELMM model with the spectral variability given by a diagonal matrix of scaling factors. The regularization imposes that, for each pixel in the image, the pixels in its neighbourhood should have the same materials and in the same proportions, up to a modulation of their spectra given by different local scaling factors. This suggests setting B to a small value, since the regularization should be valid only for small neighbourhoods. In practice, these constraints impose smoothness on the spatial factors. Therefore, the scaling factors defined by the factor matrix C could be understood as local patterns of spectral variabilities affecting each material in the image.

EXPERIMENTAL VALIDATION

Synthetic images

We chose two sets, of three reference endmembers each, corresponding to signatures of minerals from the United States Geological Survey (USGS) spectral library, comprising 224 spectral bands in the visible and near-IR. The first dataset, S 1 , presents no correlated endmembers, while the second dataset, S 2 , has two endmembers that are highly correlated. For each dataset, we simulated 10 synthetic hyperspectral images as follows. We generated 200 × 200 abundance maps using Gaussian Fields complying with the abundances sum-to-one constraint. Note that these abundance maps comprise only one pure pixel for each material, and around 5% of the pixels have an abundance coefficient superior to 0.9 for any material. We also generated spectral variability maps for each endmember using mixtures of Gaussians. Then, the pixel-dependent endmember instances were generated by multiplying the references by the corresponding spectral variability scaling factors (the achievable values are chosen so that no reflectance value becomes higher than 1, and so that the scaling factors range from 0.75 to 1.25). Then for each pixel, the mixture was performed using the linear mixing model [START_REF] Goetz | Imaging spectrometry for earth remote sensing[END_REF], and finally we added a 30dB white Gaussian noise to the image. The patch-tensor is built using a 5 × 5 sliding window, that is, B = 25.

Methodology and results

For each dataset and synthetic image, we run 20 Monte Carlo runs of the ProCo-ALS algorithm on the patch tensor, ensuring random initial factors hold the nonnegativity and, for the spatial factor, the abundances sum-to-one constraint. For the sake of comparison, we also run 20 Monte Carlo runs of the Vertex Component Analysis [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF] endmember extraction algorithm to estimate the spectral signatures from the original image, keeping the ones with highest volume. Then, we solved the FCLSU (2) to estimate the abundances. We compared the results to the actual endmembers and abundances, using the angular distance and the average root mean squared error, respectively. In Fig. 2, we compared the worst case ProCo-ALS and VCA/FCLSU results. The VCA is not very affected by the spectral variability, and since there are pure pixels in the image, the estimated endmembers are very close to the actual ones. However, the abundance estimates obtained by FCLSU could be seriously affected yielding to poor results. The ProCo-ALS obtains better estimates of the abundances, with a small loss in accuracy in the endmembers estimation.

The effect of the endmembers correlation is clear when looking at the results for the two datasets, S 1 and S 2 . On one hand, the presence of correlated endmembers yields reduced dependency on the initialization conditions. On the other hand, abundance estimates are better when data have no correlated endmembers.

Fig. 3 shows a comparison of the best results obtained by the ProCo-ALS among the 20 random initializations to the VCA/FCLSU, in terms of either endmembers or abundances estimation accuracy for the uncorrelated dataset, S 1 , and the correlated one, S 2 . For the former, the estimation of the abundances is always good, and the estimation of the endmembers is good enough (less than 5 degrees). For the latter, the accuracy of abundances depends on the accuracy of endmember estimation, and in some cases, the high correlation of the actual endmembers makes their accurate estimation a too difficult problem for a completely blind approach.

Discussion

Preliminary results show the potential of the proposed CP decomposition of hyperspectral patch-tensors for solving a blind spectral unmixing, in the presence of spectral variability. Still, some issues should be further addressed. The presence of correlated endmembers reduces the capacity of the CP decomposition to correctly identify the materials and estimate their fractional abundances. Also, it is not clear how to choose the best solution among many random initializations in practical situations, since the tensor reconstruction error appears to be uncorrelated to endmember and abundance estimation accuracy. Moreover, we are working on a modification of the proposed method to relax the regularization terms in (16), and to propose a semi-blind CP decomposition, where the spectral factors are regularized to improve their identifiability.

CONCLUSIONS

The present paper introduces two contributions. On one hand, we have proposed a patch-tensor representation of hyperspectral images, by stacking patches of the image in the third way of the tensor. On the other hand, We showed the equivalence between the CP decomposition of the proposed hyperspectral patch-tensor to a regularized version of the extended linear mixing model (ELMM) for spectral unmixing. Thus, the proposed approach aims at solving a blind spectral unmixing problem in the presence of spectral variability. We provided some preliminary results using synthetic data, and highlighted some issues and possible solutions. Further work will focus on relaxing the regularization terms of the CP decomposition equivalent ELMM, a semi-blind approach, and on providing a criterion to choose a good initialization. Also, we are working on giving evidence using realistic synthetic images and real datasets.

Fig. 2 .

 2 Fig. 2. . Box plot of the 20 runs worst case results of the CP decomposition using ProCo-ALS and random initializations, compared to the results obtained by VCA/FCLSU (in black): (a) Spectral angular distance for S 1 , (b) Abundances average RMSE for S 1 , (c) Spectral angular distance for S 2 , and (b) Abundances average RMSE for S 2 .

Fig. 3 .

 3 Fig. 3. Abundances accuracy versus endmembers accuracy for the best CP results among the 20 random initializations, using ProCo-ALS, in terms of abundances accuracy (red squares) and endmembers accuracy (green diamonds), compared to the VCA/FCLSU (blue dots): (a) Uncorrelated dataset S 1 , (b) Correlated dataset S 2 .

This work has received funding from the European Research Council under the European Community's Seventh Framework Programme FP7/2007-2013 Grant Agreement no. 320594, "DECODA".