Na Li
email: li.na@irisa.fr

Arnaud Martin
email: arnaud.martin@irisa.fr

Rémi Estival
email: remi.estival@total.com

Combination of supervised learning and unsupervised learning based on object association for land cover classification

Keywords: Object association, land cover classification, combination of supervised learning and unsupervised learning

Combination of supervised learning and unsupervised learning based on object association for land cover classification

Na Li, Arnaud Martin, Rémi

I. INTRODUCTION

As an essential mid-step for environment monitoring and change detection, the generation of thematic land cover maps through remote sensing data has been considered as one of the most important applications of remote sensing. Many studies on classification methods have been carried out to produce a reliable thematic land cover map.

Supervised learning provides thematic land cover maps with high accuracy thanks to the prior knowledge used in the training process. The prior information utilized in supervised learning can be either derieved in pixel level [START_REF] Deilmai | Comparison of two classification methods (mlc and svm) to extract land use and land cover in johor malaysia[END_REF], [START_REF] Taufik | Land cover classification of landsat 8 satellite data based on fuzzy logic approach[END_REF] or object level [START_REF] Ma | Training set size, scale, and features in geographic object-based image analysis of very high resolution unmanned aerial vehicle imagery[END_REF], [START_REF] Ma | Evaluation of feature selection methods for object-based land cover mapping of unmanned aerial vehicle imagery using random forest and support vector machine classifiers[END_REF], [START_REF] Pea-Barragn | Object-based crop identification using multiple vegetation indices, textural features and crop phenology[END_REF]. Numerous supervised classification algorithms have been developed, many of which yield very promising results. Maximum Likelihood Classification (MLC) is the most commonly used supervised classification, which assumes that training data are normally distributed. It renders satisfying performance in many studies [START_REF] Sun | Automatic remotely sensed image classification in a grid environment based on the maximum likelihood method[END_REF], [START_REF] Sisodia | Analysis of supervised maximum likelihood classification for remote sensing image[END_REF]. MLC has been proved robust for remote sensing images, as long as the data meets the distribution assumption (e.g. a Gaussian distribution) [START_REF] Wang | Comparison of ikonos and quickbird images for mapping mangrove species on the caribbean coast of panama[END_REF]. Even though an appropriate classification methodology and proper data sets are employed to generate land cover maps, supervised approaches usually require a large amount of high quality reference data for training process. However, it is not always possible to collect ground truth especially for a large or inaccessible area.

Unsupervised learning is another important group of methods applied on land cover classification [START_REF] Wang | An unsupervised classification method of remote sensing images based on ant colony optimization algorithm[END_REF]. The most highlight advantage is that it can automatically convert raw image data into distinguished clusters without prior ground knowledge. Kmeans and Iterative Self-Organizing Data Analysis Technique Algorithm (ISODATA) are the most commonly used unsupervised approaches in remote sensing thanks to their promisng performance in many studies [START_REF] Abbas | K-means and isodata clustering algorithms for landcover classification using remote sensing[END_REF], [START_REF] Li | Parallel isodata clustering of remote sensing images based on mapreduce[END_REF]. However, even though unsupervised approaches are more appropriate for large and inaccessiable areas, they can not directly produce thematic land cover maps; therefore time-consuming manually labeling process is always required. Due to the lack of prior ground information, unsupervised learning is usually less accurate than supervised learning [START_REF] Ismail | Evaluating supervised and unsupervised techniques for land cover mapping using remote sensing data[END_REF].

As both supervised and unsueprvised methods have their pros and cons, it is thus necessary to consider to combine them together for better results. Although many studies have been carried out to achieve the combination, most of them focus on hierarchical structures which makes supervised learning and unsupervised learning dependent [START_REF] Markowska-Kaczmar | Combined unsupervisedsupervised classification method, in Knowledge-Based and Intelligent Information and Engineering Systems[END_REF]. This kind of combination is less robust as the mistakes from the former classifier may be hard to reduce by the latter especially when it is unsupervised. Many combinations have been also developed in feature level [START_REF] Castelli | Combination of supervised and unsupervised learning for training the activation functions of neural networks[END_REF] to generate more representative feature spaces or hyper planes; therefore it is difficult to measure uncertainties of the results from each individual classifier.

The combination in the decision level of independent supervised and unsupervised classifiers is still a big challenge. In [START_REF] Karem | Combination of supervised and unsupervised classification using the theory of belief functions[END_REF], the authors applied Dempster-Shafer theory to achieve this combination, considering supervised learning and unsupervised learning as independent sources, through which the uncertainty and the imprecision in the fusion can be measured. However, currently existing methods do not focus on the conflict situation where the number of supervised classes is not the same as the number of unsupervised clusters.

Therefore, in this paper, we proposed a novel framework using object association based on Dempster-Shafer theory to This paper is organized as follows: In section II, we introduce the study area and data set, followed by some bases of Dempster-Shafer theory in section III. Object association problem is explained in section IV. Section V presents the proposed framework in details. Then, section VI illustrates results of our experiments. Conclusions are drawn in section VII.

II. STUDY AREA AND DATA SETS

A. Study area

The area of study located in Colorado in United Status, lies within latitudes 38 • 30 53.44 N -36 • 22 23.05 N and longitudes 108 • 0 10.24 W -105 • 18 35.60 W. It is contained San Juan National Forest and Rio Grande National Forest, as shown in Fig. 1. The two national forest parks are mainly composed of deciduous forest, evergreen forest, mix forest. The remaining natural vegetation still occurs as shrub, grassland and herbaceous. A large area of pasture is also contained near Monte Vista. This selected area contains two national forest parks as well as various developed areas thus can be used to mimic the complex environment with a wide range of land covers. We separate the study area into two subareas of equal size: validation area and training area.

B. Satellite

The satellite used in this study is LandSat-8 OLI consisting of eight spectral bands with a spatial resolution of 30 meters, a panchromatic band with a resolution of 15 meters and two thermal bands with a resolution of 100 meters. The image acquired on the 11 June 2018 was obtained from USGS Earth Explorer [16]. Geometric correction of the image has already done through UTM map projection by NASA. The bands used for our study include band 1-7 which are all spectral bands.

C. Training and validation data

Training and validation data applied in our study was obtained from National Land Cover Database 2011 (NLCD 2011) [START_REF] Homer | Completion of the 2011 national land cover database for the conterminous united statesrepresenting a decade of land cover change information[END_REF], which was created by the Multi-Resolution Land Characteristics (MRLC) Consortium. NLCD 2011 has 16class land cover classification scheme that has been applied consistently across the United States at a spatial resolution of 30 meters. The code of the land cover classification scheme of NLCD 2011 is shown in Table I. In our study, we select half of the study area as the training area which includes all 16 types of land covers mentioned in Table I with the size of 2000×4000 pixels to generate training samples, as shown in Fig. 1. The other half of the study area with the same size is chosen as validation area as shown in Fig. 1. To simulate the situation where training samples are insufficient for supervised classifier, only a small amount of pixels per class are randomly selected as training samples, which will be explained in details in section VI

III. DEMPSTER-SHAFER THEORY

As a generalization of traditional probability, Dempster-Shafer theory [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF], [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF] allows to distribute support to not only a single proposition itself but also to the union of propositions. One of the greatest advantages of Dempster-Shafer theory is that it allows to take into consideration uncertainty and imprecision at the same time via two functions: belief and plausibility, derived from mass function. The mass function is defined on all the subsets of the frame of discernment Ω = {ω 1 , . . . , ω n }, and assigns belief degree to all the elements in the power set of discernment, noted as 2 Ω .

The sum of the masses of all the propositions is one:

A⊆Ω m(A) = 1 (1)
On the contrary, the mass function of the null proposition ∅ is usually set to zero but it is also possible to be a positive value. m(∅) = 0 refers to a closed world hypothesis, in which the discernment Ω contains all the possible situations in reality.

If m(∅) is superior to 0, this corresponds to an open world hypothesis, where unknown situation outside of Ω can be considered. m(A) can be considered as a degree of evidence supporting the claim that a specific element of Ω belongs to the set A, yet not to any subset of A [START_REF] Klir | Fuzzy Sets, Uncertainty, and Information[END_REF].

IV. OBJECT ASSOCIATION

Object association was originally derived from tracking problem to estimate the status of mobile objects detected by multiple sensors. Mathematically speaking, object association refers to match two finite set of objects E = {e 1 , ..., e n } and F = {f 1 , ...f p } with possibly different combinations. It is usually assumed that each object in one set should be matched with at most one object in the other set. An object is possible to disappear between two successive time frames, which leads to no correspondent counterpart for this object in the other set.

The problem could be considered in the framework of Dempster-Shafer theory, also referred as Belief function. Therefore, a piece of evidence presenting the possible association of two objects e i and f j can be modeled by a mass function m ij on the discernment Θ ij = {0, 1}. Let R ij represent the relation between the object e i and f j . The mass function m ij ({1}) refers to the probability of R ij = 1, m ij ({0}) is the probability R ij = 0 and m ij ({0, 1}) is the probability of knowing nothing about the relation R ij .

Several solutions about choosing the best relation R * based on pairwise mass functions m ij have been proposed. In [START_REF] Mercier | Object association with belief functions, an application with vehicles[END_REF], the authors combines the mass function m ij for each i and then find the best relation R * through maximum pignistic probability. In [START_REF] Denoeux | Optimal object association in the dempster-shafer framework[END_REF], the problem is transferred into an equivalent linear assignment problem and then searches for the most plausible relation as the best one, which thus can be solved in polynomial time. Since the number of classes and clusters is large in our study, we applied the approach proposed in [START_REF] Denoeux | Optimal object association in the dempster-shafer framework[END_REF] on the combination of supervised learning and unsupervised leaning. It can be stated as the following inter linear program:

max( n i=1 p j=1 ω ij R ij ) (2)
subject to

p j=1 R ij ≤ 1 ∀i ∈ {1, ..., n} (3)
n i=1 R ij ≤ 1 ∀j ∈ {1, ..., p} (4)
R ij ∈ {0, 1} ∀i ∈ {1, ..., n}, ∀j ∈ {1, ..., p} (5)
in which w ij is a coefficient generated by m ij presenting the uncertainty as well as imprecision on the evidence of relation between e i and f j .

V. THE PROPOSED FRAMEWORK A. Combination of supervised learning and unsupervised learning Let us denote S as the set of classes generated from supervised learning and C as the set of clusters derived from unsupervised learning. Based on the principle of object association, we have to match elements from two sets of objects S = {s 1 , ..., s n } and C = {c 1 , ..., c p } with the most possible combination.

As explained in Section IV, mass function m ij with 1 ≤ i ≤ n, 1 ≤ j ≤ q, indicates all available evidences about the association between S and C. Each m ij encodes a piece of evidence about a binary variableR ij . This variable indicates s i is the same type as c j when it equals to 1, and 0 otherwise. Each m ij is thus on the discernment Θ = {0, 1} having:

m ij ({1}) = α ij (6) m ij ({0}) = β ij (7)
m ij ({0, 1}) = 1 -α ij -β ij (8)
where α ij and β ij is measured based on the similarity between a class s i and a cluster c j . We assumed that n < p because the number of class, a kind of prior knowledge, can be consider as the lower bound of the number of cluster which is an unknown variable in real situation. Hence, a cluster c j is usually smaller than a class s i . To measure the similarity between s i and c j , we consider the percentage of pixels in c j is labeled as s i . We thus use the ratio of pixels in intersection of s i in c j to pixels in c j as the similarity δ ij , shown as follows:

δ ij = |s i ∪ c j | |c j | (9)
For α ij and β ij , we have:

α ij = γδ ij ( 10
)
β ij = γ(1 -δ ij ) (11)
in which γ represents a discounting coefficient ranging from 0 to 1 to measure the reliability of classifiers.

To assign label of classes to clusters, we consider the following linear optimization problem:

max R ( i,j ω ij R ij ) (12)
Fig. 2: The proposed framework.

subject to

p j=1 R ij ≥ 1 ∀i ∈ {1, ..., n} (13)
n i=1 R ij = 1 ∀j ∈ {1, ..., p} (14)
R ij ∈ {0, 1} ∀i ∈ {1, ..., n}, ∀j ∈ {1, ..., p} (15)
with

ω ij = ln 1 -β ij 1 -α ij (16)
The most plausible relation R * , which indicates assignments between classes and clusters, can thus be found by solving the above optimization problem. Unlike to classical object association problem where each object in one set has at most one counterpart in the other set, in our problem, each class s i can be assigned to multiple clusters while each cluster c j is merely matched with one class as we always have n < p.

B. Generation of land cover map for a large area

The proposed framework can produce a thematic land cover map on a large area merely by deploying a supervised approach on a small slice of the validation region. We select a slice of area containing all types of land covers to employ a supervised learning approach. On the contrary, unsupervised learning approach is applied on the whole area. The principal ideal behind the proposed framework is to assign the most plausible label to each clusters in the same area where thematic labels are offered by the supervised approach; later we can enlarge these labeled clusters to the whole area.

The main steps of the proposed framework is shown in Fig 2, in which the region marked with orange rectangular denoted as A is the selected slice for conducting supervised learning approach. The set of classes produced by the supervised approach deployed on the region A is denoted as:

S A = {s A 1 , ..., s A n } (17)
The set of clusters from the unsupervised learning on the entire validation area is denoted as:

C = {c 1 , ..., c p } (18)
The set of clusters in the region A denoted as C A is defined as:

C A = {c A i |c A i = c i ∩ S A } (19)
The object association is firstly applied between the set S A and the set C A ; therefore clusters in C A can be all assigned by a thematic label from the supervised approach. To enlarge thematic labels to the whole area representing by the set of clusters C, we assign each cluster c i with the same label as its subset c A i shown as follows:

L(c i ) = L(c A i ) (20)
where L(x) represents the final thematic label of object x.

VI. EXPERIMENT RESULTS AND DISCUSSION

In this section, we tested the proposed framework with different combinations of supervised learning and unsupervised learning. To verify the possibility to apply the proposed framework on a large area, we only carried out supervised learning on a quarter of the validation area. Based on the labels it generates, a thematic land cover map of the entire validation area can be produced by the proposed framework, as illustrated in Fig. 2. We also applied supervised learning on the entire validation area to compared with the results from the proposed framework.

Apart from the experiments on different combination of methods, we also conducted the proposed framework with different size of training data and made the comparison with supervised learning. We repeated all the experiments in our study 10 times to take their average results.

A. Combination of MLC and ISODATA

We have first used the combination of MLC and ISODATA in the proposed framework with the number of clusters ranging from 20 to 500. To mimic the situation where training samples are severely insufficient, we first selected a very small amount of training samples to train MLC. Table II summarizes the comparison between the proposed framework and MLC on the validation area with training samples accounting for 0.2% in the training area. The overall accuracy and the accuracy of each land cover category in Table II were measured by confusion matrix. Due to the limitation of space, we have employed the code of NLCD classification mentioned in Table I to represent each land cover category and marked the best accuracy in bold. The overall accuracy of the proposed framework improves as the number of clusters k increasing and the best performance is reached when k = 100.

We displayed the results of MLC shown in Fig. 3a to compare with the best thematic land cover map generated through the combination of MLC and ISODATA, as shown in Fig. 3b. MLC apparently misclassifed a part of shrub as medium intensity developed area. However, after the combination with ISODATA in the proposed framework, the best Fig. 3: Thematic land cover maps accuracy of shrub is enhanced to 0.684%. This improvement is presented by the reduction of the red area in Fig. 3b.

Besides shrub, other terrestrial vegetation such as deciduous forest, evergreen forest and grassland also have good improvements of classification accuracy after the combination. However, some land cover categories such as developed areas, barren land and emergent herbaceous wetlands, shows unsatisfying classification accuracy, which are also difficult to We also tested the proposed framework and MLC with the different quantities of training samples fixing k = 100. We started the test on the severely insufficient training samples which only represent 0.05% of the training area. Table . III reports the overall accuracy of MLC and the proposed framework with the combination of MLC and ISODATA on the validation area. Fig. 4. gives a more apparent display of the enhancement of their accuracies with the increasing of training samples. The proposed framework always performs better than MLC when the percentage of training samples ranging from 0.05% to 100% and the largest enhancement of accuracy occurs at 0.2%.

B. Combination of SVM and ISODATA

A combination of Support Vector Machine (SVM) and ISODATA was also considered in the proposed framework. In our study, we selected Radial Basis Function (RBF) kernel as it is proved to have better results for LandSat data with 7 bands [START_REF] Huang | An assessment of support vector machines for land cover classification[END_REF]. The error weight and Gamma value in SVM were set separately to 100 and 0.167 as is proposed in [START_REF] Yang | Parameterizing support vector machines for land cover classification[END_REF].

Table IV reports the overall accuracy as well as the accuracy of each individual land cover category when training samples account for 0.2%. The performance of SVM on the entire validation area is also displayed at the bottom of the Table IV to make the comparison. The number of clusters k ranged from 20 to 500 with the best performance of each land cover category marked in bold. For this combination, the best overall accuracy is reach when k is 60 and the correspondent results are shown in Fig. 5a and Fig. 5b.

Compared with MLC, SVM is obviously less competitive to distinguish shrub and medium intensity developed area when training data is severely insufficient. This kind of combination has a good improvement of the accuracy of shrub and grassland. However, the enhancement of overall accuracy from 0.441 to 0.467 is less apparent than the combination of MLC and ISODATA.

Another experiment to exterminate the performance of the proposed framework with different size of training data was also conducted. The number of clusters k in ISODATA is fixed at 60 at which this combination achieves the best performance in the previous experiment. The results are summarized in Table V and Fig. 6. The overall accuracy is improved when percentage of training data is from 0.05% to 0.8%, which is less satisfying than the combination of MLC and ISODATA.

C. Discussion

We applied two combinations in the proposed framework: MLC with ISODATA, and SVM with ISODATA. For the experiments with different number of clusters k, we all started when k is 20, which is approximate to the number of classes set by MLC and SVM. However, the two combinations are both less accurate than MLC or SVM in this situation because of the limitation of unsupervised learning. As a clustering method, ISODATA separates data merely based on their spectral dissimilarities in the seven optical bands used in the study. In the complex natural environment, however, different types of land cover may have similar spectral properties and the same type may reflects differently in spectral bands due to varieties of temperature and illumination. A direct application of ISODATA on our study data is thus less satisfying. The proposed framework based on object association is to achieve the combination in the decision level. It can be considered as an average results of supervised learning and unsupervised learning when the number of clusters is approximate to the number of classes. Therefore, the combination of MLC and ISODATA or SVM and ISODATA is less accurate than the original supervised approach. A relatively large value of k in ISODATA can separate original data into multiple homogeneous clusters which thus can be considered as over-segmentation or super pixels in spectral space. In this situation, pixels in each pure cluster are gathered as an entire object whose label is associated by the proposed framework. Compared to an individual supervised approach in pixel level, the proposed framework can generate the thematic land cover maps on object level; therefore, it has an evident improvement of the overall accuracy. With the increment of k, the combination of SVM and ISODATA is less efficient than MLC and ISODATA probably because the spectral properties of land covers in our study area are more distinguishable for MLC than SVM when training data is insufficient.

With the different size of training data, the combination of MLC and ISODATA always renders better performances than MLC alone and this improvement of accuracy is most satisfying when the training samples account for 0.2% in the training area. However, for SVM and ISODATA, the enhancement of overall accuracy only occurs when training samples are from 0.1% to 0.8 possibly because SVM is hard to find appropriate hyper planes for classification through such limited training data. Therefore, this combination is less effective than the first one in our study area.

The uncertainty and imprecision are measured through mass functions in the object association. We applied the ratio of pixels in intersection of a class s i in a cluster c j to pixels in c j as the similarity between a class and a cluster. The uncertainty thus reduces when the number of clusters increases, which explains why results become better with large k. However, when k becomes too large, the results are less competitive due to the loss of information in unsupervised learning.

VII. CONCLUSION

We proposed a novel framework to combine supervised learning and unsupervised especially focusing on the conflict situation where the number of classes is not the same as the number of clusters. It can thus be applied in large or inaccessible areas where not enough ground truth is available.

The combination conducted in the framework takes merely spectral information as consideration. Therefore, in our future work, we will concentrate on enlarging this framework to combine various formats of information such as contextual information and ancillary data in the decision level to generated more accurate thematic land cover maps. On the other hand, as uncertainty and imprecision take important roles in classification, the application of Dempster-Shafer theory in the fusion framework is also another promising aspect in our future work.

Fig. 1 :

 1 Fig. 1: Map of the study area

  (a) Results of MLC. (b) Results of MLC + ISODATA when k=100.

Fig. 4 :

 4 Fig. 4: Accuracy with different size of training samples.

  (a) Results of SVM. (b) Results of SVM + ISODATA with k=60.

Fig. 5 :

 5 Fig. 5: Thematic land cover maps

Fig. 6 :

 6 Fig. 6: Accuracy with different size of training samples

TABLE I :

 I Land cover classification scheme of NLCD 2011

	Code Legend	Code Legend
	11	Open Water	42	Evergreen Forest
	12	Perennial Ice/Snow	43	Mixed Forest
	21	Developed, Open Space	52	Shrub/Scrub
	22	Developed, Low Intensity	71	Grassland/Herbaceous
	23	Developed, Medium Intensity	81	Pasture/Hay
	24	Developed, High Intensity	82	Clutivated Crops
	31	Barren Land (Rock/Sand/Clay)	90	Woody Wetlands
	41	Deciduous Forest	95	Emergent Herbaceous Wetlands

TABLE II :

 II Comparison of MLC and MLC + ISODATA in classification accuracy.

	Number of						Accuracy of each land cover category								Overall
	clusters	11	12	21	22	23	24	31	41	42	43	52	71	81	82	90	95	accuracy
	20	0.316 0.051 0.078 0.077 0.133 0.019 0.064 0.212 0.494 0.464 0.352 0.205 0.255 0.142 0.116 0.199 0.325
	40	0.590 0.031 0.041 0.057 0.064 0.008 0.034 0.189 0.672 0.574 0.571 0.452 0.257 0.163 0.047 0.235 0.475
	60	0.541 0.001 0.026 0.079 0.044 0.034 0.037 0.179 0.694 0.572 0.634 0.451 0.306 0.075 0.039 0.267 0.496
	80	0.463 0.043 0.027 0.063 0.025 0.022 0.059 0.198 0.711 0.521 0.684 0.429 0.308 0.071 0.027 0.294 0.505
	100	0.579 0.026 0.019 0.069 0.033 0.032 0.046 0.206 0.698 0.606 0.676 0.479 0.313 0.048 0.014 0.339 0.514
	200	0.746 0.006 0.027 0.088 0.098 0.012 0.079 0.164 0.682 0.595 0.659 0.453 0.319 0.065 0.040 0.347 0.501
	300	0.697 0.019 0.020 0.088 0.077 0.022 0.0941 0.232 0.684 0.580 0.672 0.441 0.318 0.057 0.046 0.363 0.504
	400	0.689 0.005 0.021 0.092 0.115 0.0115 0.084 0.219 0.678 0.600 0.657 0.441 0.324 0.091 0.047 0.352 0.500
	500	0.739 0.012 0.035 0.110 0.106 0.057 0.101 0.239 0.664 0.551 0.624 0.469 0.307 0.119 0.062 0.367 0.497
	MLC	0.771 0.064 0.065 0.109 0.375 0.294 0.591 0.231 0.572 0.764 0.449 0.377 0.278 0.262 0.180 0.327 0.432

TABLE III :

 III Overall accuracy of MLC and MLC+ISODATA with different size of training samples.

	Percentage
	in the training	0.05% 0.1%	0.2%	0.4%	0.8%	1.2%	1.6%	2.0%	2.4%	5.0%	10.0% 20.0% 40.0% 80.0% 100.0%
	area
	MLC	0.418	0.437	0.446	0.447	0.451	0.453	0.456	0.457	0.456	0.460	0.470	0.484	0.499	0.528	0.539
	MLC+ISODATA	0.475	0.504	0.513	0.510	0.505	0.508	0.507	0.509	0.509	0.511	0.514	0.518	0.529	0.554	0.562

TABLE IV :

 IV Comparison of SVM and SVM + ISODATA in classification accuracy

	Number of						Accuracy of each land cover category								Overall
	clusters	11	12	21	22	23	24	31	41	42	43	52	71	81	82	90	95	accuracy
	20	0.010 0.070 0.202 0.104 0.245 0.001 0.219 0.285 0.569 0.436 0.193 0.126 0.263 0.161 0.202 0.232 0.321
	40	0.014 0.011 0.056 0.046 0.562 0.007 0.274 0.243 0.661 0.270 0.126 0.584 0.268 0.149 0.112 0.084 0.445
	60	0.737 0.032 0.006 0.097 0.437 0.001 0.320 0.411 0.629 0.003 0.503 0.424 0.344 0.029 0.045 0.214 0.467
	80	0.756 0.187 0.026 0.092 0.491 0.006 0.351 0.406 0.668 0.373 0.335 0.427 0.278 0.073 0.058 0.328 0.457
	100	0.772 0.027 0.041 0.124 0.565 0.008 0.383 0.336 0.670 0.227 0.223 0.464 0.312 0.053 0.092 0.287 0.450
	200	0.739 0.025 0.015 0.153 0.545 0.001 0.339 0.342 0.665 0.262 0.214 0.447 0.347 0.014 0.149 0.299 0.448
	300	0.739 0.008 0.015 0.153 0.545 0.001 0.339 0.342 0.665 0.262 0.214 0.447 0.346 0.014 0.149 0.299 0.443
	400	0.681 0.010 0.033 0.169 0.649 0.002 0.368 0.331 0.684 0.401 0.232 0.422 0.331 0.113 0.125 0.312 0.450
	500	0.716 0.014 0.013 0.207 0.541 0.001 0.410 0.318 0.668 0.404 0.258 0.423 0.338 0.094 0.146 0.335 0.451
	SVM	0.765 0.058 0.057 0.202 0.686 0.106 0.503 0.313 0.683 0.522 0.202 0.439 0.321 0.209 0.158 0.359 0.441

TABLE V :

 V Overall accuracy of SVM and SVM+ISODATA with different size of training samples.

	Percentage
	in the training	0.05% 0.1%	0.2%	0.4%	0.8%	1.2%	1.6%	2.0%	2.4%	5.0%	10.0% 20.0% 40.0% 80.0% 100.0%
	area
	SVM	0.381	0.428	0.442	0.461	0.464	0.490	0.493	0.493	0.492	0.525	0.571	0.604	0.676	0.709	0.713
	SVM+ISODATA	0.392	0.423	0.468	0.470	0.472	0.485	0.483	0.477	0.482	0.516	0.556	0.586	0.658	0.681	0.679


Hsin-Ying Lee
email: hylee@ee.ncku.edu.tw

Yu-Chang Lin

Yu-Ting Su

Chia-Hsin Chao

Véronique Bardinal

Prof Ying Hsin-

Lee

Performance improvement of GaN-based flip-chip white light-emitting diodes with diffused nanorod reflector and with ZnO nanorod antireflection layer

Keywords: Aqueous solution method, flip-chip white light-emitting diodes, diffused nanorod reflector, ZnO nanorod

The GaN-based flip-chip white light-emitting diodes (FCWLEDs) with diffused ZnO nanorod reflector and with ZnO nanorod antireflection layer were fabricated.

The ZnO nanorod array grown using an aqueous solution method was combined with Al metal to form the diffused ZnO nanorod reflector. It could avoid the blue light emitted out from the Mg-doped GaN layer of the FCWLEDs, which caused more blue light emitted out from the sapphire substrate to pump the phosphor. Moreover, the ZnO nanorod array was utilized as the antireflection layer of the FCWLEDs to reduce the total reflection loss. The light output power and the phosphor conversion efficiency of the FCWLEDs with diffused nanorod reflector and 250-nm-long ZnO nanorod antireflection layer were improved from 23.90 mW to 21.15 mW and 77.6% to 80.1% in comparison with the FCWLEDs with diffused nanorod reflector and without ZnO anaorod antireflection layer, respectively.

Introduction

Recently, people have put more and more attentions on the problems of energy shortage and environment problems. Many research teams actively investigated how to improve energy saving and how to product renewable energy. The white light-emitting diodes (WLEDs) are expected to be the promising green lighting sources and have been widely used in solid-state lighting and taken as the next generation lighting sources to replace the compact fluorescent lamp (CFL) and incandescent light. The WLEDs exhibit many superior advantages such as the long lifetime, high luminous intensity, energy saving, fast response, low heat dissipation and high reliability. In general, there are three approaches can be used to fabricate WLEDs: (1) using the blue light-emitting diodes (LEDs) to excite the yellow phosphor [START_REF] Rao | An improved slurry method of self-adaptive phosphor coating for white pc-LED packaging[END_REF][START_REF] Lee | Light distribution and light extraction improvement mechanisms of remote GaN-based white light-emitting-diodes using ZnO nanorod array[END_REF], (2) using an ultraviolet (UV) LEDs to excite red, green, and blue phosphors [START_REF] Mao | Tricolor emission Ca 3 Si 2 O 7 :Ln (Ln=Ce, Tb, Eu) phosphors for near-UV white light-emitting-diode[END_REF][START_REF] Park | Melilite-structure CaYAl 3 O 7 :Eu 3+ phosphor: structural and optical characteristics for near-UV LED based white light[END_REF], and (3) mixing three primary colors of red, green, and blue LEDs [START_REF] Muthu | Red, green, and blue LEDs for white light illumination[END_REF]. However, the large difference of refractive index between the GaN-based semiconductor (n GaN =2.4) and the air (n air =1) led to the increase of total internal reflection and the Fresnel loss. Generally, several methods, such as die shape [START_REF] Kuo | Nitride-based LEDs with phosphoric acid etched undercut sidewalls[END_REF][START_REF] Brinkley | Chip shaping for light extraction enhancement of bulk c-plane light-emitting diodes[END_REF],

surface roughness [START_REF] Lee | Light output enhancement of GaN-based roughened LEDs using bias-assisted photoelectrochemical etching method[END_REF][START_REF] Lin | Enhanced light extraction of GaN-based light emitting diodes using nanorod arrays[END_REF], photonic crystal [START_REF] Kim | Enhanced light extraction from GaN-based light-emitting diodes with holographically generated two-dimensional photonic crystal patterns[END_REF][START_REF] Cho | Light extraction enhancement from nano-imprinted photonic crystal GaN-based blue light-emitting diodes[END_REF], reflector [START_REF] Chang | Improved light extraction efficiency of a high-power GaN-Based light-emitting diode with a three-dimensional-photonic crystal (3-D-PhC) backside reflector[END_REF][START_REF] Lobo | Enhancement of light extraction in ultraviolet light-emitting diodes using nanopixel contact design with Al reflector[END_REF], and antireflection coating layer [START_REF] Yan | Deep-ultraviolet tailored-and low-refractive index antireflection coatings for light-extraction enhancement of light-emitting diodes[END_REF], have been reported previously to overcome the undesired total internal reflection and Fresnel loss. Among these methods, the reflector exhibits several advantages including simple fabrication, low cost, and outstanding light extraction improvement of LEDs. In this work, to improve the light extraction efficiency of WLEDs, the designed diffused nanorod reflector constructed with the 700-nm-thick Al metal and ZnO nanorod array with optimal rod-length of 500 nm was used in the flip-chip white light-emitting diodes (FCWLEDs). The detail investigation of diffused nanorod reflector was shown in previous published paper [START_REF] Lee | Light extraction enhancement of GaN-based flip-chip light-emitting diodes using diffused nanorod reflector[END_REF]. The diffused nanorod reflector could effectively increase the reflection angle and enhance the light escaping probability from the semiconductor compared with the traditional flat reflector. Furthermore, to further improve the light extraction from the light-emitted side of the FCWLEDs, the ZnO nanorod arrays were also used as the antireflection layer. Finally, the yellow phosphor layer was coated on the ZnO nanorod antireflection layer by using the remote phosphor coating technique which possessed the uniformity and high color converted efficiency.

Experiment Procedure

Figure 1 shows the schematic configuration of the GaN-based FCWLEDs with diffused nanorod reflector and with ZnO nanorod antireflection layer. A metal organic chemical vapor deposition (MOCVD) system was used to grow the epitaxial layers of the GaN-based blue LEDs. A 2.8-μm-thick GaN buffer layer, a 4-μm-thick Si-doped GaN layer (n=310 17 cm -3 ), an undoped InGaN/GaN multiple quantum well (MQW) active layer, a 33-nm-thick Mg-doped Al 0.2 Ga 0.8 N layer (p=110 17 cm -3 ), and a 150-nm-thick Mg-doped GaN layer (p=310 17 cm -3 ) were sequentially grown on the c-plane sapphire substrates. The InGaN/GaN MQW active layer was consisted of twelve pairs of 3-nm-thick In 0.2 Ga 0.8 N well layer and 12.7-nm-thick GaN barrier layer.

The fabrication processes of the conventional GaN-based LEDs are as follows. The epitaxial layer with mesa (300×300 μm 2 ) pattern protected by a 300-nm-thick Ni metal mask was etched down to the Si-doped GaN layer by a reactive ion etching (RIE) system. After the mesa etching, the remaining Ni metal mask was removed using aqua regia. The n-electrode Ti/Al/Ti/Au (25/100/50/150 nm) was deposited using an electron-beam evaporator. Moreover, the n-electron was treated by a rapid thermal annealing (RTA) system with 850 o C for 2 min in a pure N 2 ambient for obtaining the good ohmic contact performance [START_REF] Lee | Long term thermal stability of Ti/Al/Pt/Au ohmic[END_REF]. Thin Ni/Au (2.5/2.5 nm) and thick Ni/Au (20/100 nm) metal contacts were deposited on Mg-doped GaN layer as were also fabricated for comparison.

Result and Discussion

To find the optimal rod-length of the ZnO nanorod in the antireflection layer, the performances of the GaN-based FCLEDs with diffused nanorod reflector and with 0-, 100-, 250-, 500-, and 750-nm-long ZnO nanorod arrays were firstly investigated.

Figure 2 shows the light output power-current (L-I) curves of the GaN-based FCLEDs with diffused nanorod reflector and with ZnO nanorod arrays of various rod-lengths measured using an integrating sphere. As shown in Figure 2, the light output power of the FCLEDs with diffused nanorod reflector and with 250-nm-long ZnO nanorod array was improved from 30.85 mW to 34.99 mW, compared with the FCLEDs with diffused nanorod reflector and without ZnO nanorod (i.e., 0-nm-long ZnO nanorod) antireflection layer at an injection current of 280 mA. However, when the rod-length of the ZnO nanorod array was longer than 250 nm, the light output power decreased owing to the absorption of the ZnO nanorod array. To confirm the optimal rod-length of the ZnO nanorod antireflection layer on the output side of the FCLEDs with diffused nanorod reflector, the standard optical software TracePro was also used to simulate the light propagated within the ZnO nanorod arrays of various rod-lengths.

Figure 3 shows the simulating results of an enhancement in the illumination efficiency for the FCLEDs with diffused nanorod reflector and with ZnO nanorod antireflection layer of various rod-lengths, in comparison with the FCLEDs with diffused ZnO nanorod reflector and without ZnO nanorod antireflection layer. The illumination efficiency of the FCLEDs with diffused nanorod reflector and with 100-, 250-, 500-, and 750-nm-long ZnO nanorod antireflection layer shown an increase of 10.47%, 11.97%, 9.32%, and 8.33%, respectively. The experimental results of the enhanced proportion of the light output power for the resulted FCLEDs are also shown in Figure 3. The experimental results and the simulation results had the similar tendency. In addition, the optimal effective refractive index of the antireflection layer can be determined by the formula of (n air  n GaN ) 1/2 , where n air of 1.00 and n GaN of 2.40 were the refractive index of air and n-GaN, respectively. The solution of the optimal effective refractive index of the antireflection layer was 1.55. In this work, the effective refractive index of the 250-nm-long ZnO nanrod antireflection layer was measured by an ellipseometer and the value was about 1.57 which was similar to the optimal effective refractive index of the antireflection layer. Therefore, the rod-length of 250 nm in the ZnO nanorod antireflection layer was the optimal length for the FCLEDs. where P p is the optical power emitted from the blue FCLEDs without phosphor layer and P phosphor is the optical power converted from the phosphor layer. Compared with the FCWLED B, the phosphor conversion efficiency of FCWLED A operated at an injection current of 20 mA was improved from 77.6% to 80.1%. The improvement was attributed to that the ZnO nanorod antireflection layer could reduce the light total reflection at the sapphire surface, which could guide more blue light to pump the phosphor layer. Generally, the phosphor conversion efficiency as a function of the injection current of the LEDs. Consequently, the EL spectra of the FCLED C and FCWLED A at various injection currents were measured and are shown in Figure 5 (a)

and (b), respectively. Using the experimental results and Eq. ( 1) to estimate the phosphor conversion efficiency of FCWLED A at various injection currents were measured and is shown in Figure 6. The phosphor conversion efficiency of FCWLED A was enhanced with an increased of the injection current, a maximum value of 92.5% was obtained as the FCWLED A operated at a current of 80 mA. This phenomenon was attributed to that the phosphor conversion efficiency reached saturation at an injection current of 80 mA. Therefore, when the injection current was more than 80 mA, the phosphor conversion efficiency decreased owing to the excess blue light directly emitted out the FCWLEDs or the excess photons multi-scattered with the phosphor and were absorbed inside the FCWLEDs. The L-I curves of FCWLED A and FWCLED B were also measured using an integrating sphere system.

The light output power of FCWLED A and FWCLED B was 23.90 mW and 21.15 mW, respectively, at an injection current of 280 mA.

Figure 7 shows the chromaticity diagram of FCWLED A and FCWLED B at an injection current of 80 mA. As shown in Figure 7, the chromaticity coordinate of the emission of the corresponding LEDs exhibited the same coordinate of (0.335, 0.333), which demonstrated that the emitted light was a white light.

Conclusions

The light output power and phosphor conversion efficiency of the FCWLEDs were successfully improved by using the diffused nanorod reflector and the ZnO nanorod antireflection layer. The diffused nanorod reflector constructed by the 500-nm-long ZnO nanorod arrays and the high reflective Al metal could effectively enhance the probability of light reflectivity from the bottom of the FCWLEDs. Moreover, according to the experimental results and the simulation results, the optimal 250-nm-long ZnO nanorods grown on the output side of the LEDs as the

  The commercial WLEDs were commonly integrated by the GaN-based blue LEDs with yellow phosphor (Y 3 Al 5 O 12 :Ce 3 + or YAG:Ce 3 + ) layer. To augment the application of WLEDs, the development of high performance WLEDs are immediately needed. Since the properties of the GaN-based blue LEDs seriously affect the performances of WLEDs, how to improve the internal quantum efficiency and external quantum efficiency of the blue LEDs are the primarily issues. In view of the significant progress of GaN-based epitaxial technology, the internal quantum efficiency of the GaN-based blue LEDs achieved the theoretical limit of 90%.

Figure 1 .

 1 Figure 1. Then a 50-nm-thick AZO seed layer and ZnO nanorod arrays were

Figure 4

 4 Figure4shows the electroluminescence (EL) spectra of the FCWLEDs with

Figure captions FIGURE 1 FIGURE 2 FIGURE 3 FIGURE 4 FIGURE 5 FIFURE 6 FIGURE 7

 captions1234567 Figure captions

  FIGURE 5

  antireflection layer could guide more blue light and could excite more phosphor to enhance the phosphor conversion efficiency. The light output power of the FCWLEDs with diffused nanorod reflector and with ZnO nanorod antireflection layer was enhanced 13.0% in comparison of the FCWLEDs with diffused nanorod reflector and without with ZnO nanorod antireflection layer at injection current of 280 mA. The phosphor conversion efficiency was also improved 2.5%. contacts to n-type GaN," Appl. Phys. Lett., vol. 76, no. 17, pp. 2364-2366, 2000.[17] C. S. Lee, Y. J. Lin, and C. T. Lee, "Investigation of oxidation mechanism for ohmic formation in Ni/Au contacts to p-type GaN layers," Appl. Phys. Lett., vol.

	79, no. 23, pp. 3815-3817, 2001.

Acknowledgement

This work was supported by the Industrial Technology Research Institute, Bureau of Energy, Ministry of Economic Affairs of Taiwan, Contract No. 102-E0605, the National Sciences Council of Taiwan under Grant NSC-101-2923-E-006-004-MY2, and the Advanced Optoelectronic Technology Center, National Cheng Kung University, Taiwan.


H Chraibi

A Dutfoy

T Galtier
email: tgaltier@gmail.com

J Garnier

OPTIMAL INPUT POTENTIAL FUNCTIONS IN THE INTERACTING PARTICLE SYSTEM METHOD

Keywords: Rare event, Reliability assessment, Feymann-Kac particle filters, interacting particles system, Sequential Monte-Carlo Samplers, selection, potential function 2010 Mathematics Subject Classification: Primary 60K10 Secondary 90B25, 62N05

The assessment of the probability of a rare event with a naive Monte-Carlo method is computationally intensive, so faster estimation methods, such as variance reduction methods, are needed. We focus on one of these methods which is the interacting particle (IPS) system method.

The method requires to specify a set of potential functions. The choice of these functions is crucial, because it determines the magnitude of the variance reduction. So far, little information was available on how to choose the potential functions. To remedy this, we provide the expression of the optimal potential functions minimizing the asymptotic variance of the estimator of the IPS method.

Introduction

Let (Z k ) k∈{0,...,n} be a Markov chain with values in the measurable spaces (E k , E k ), and with initial law ν 0 , and with a kernel Q k such that for k > 0 and for any bounded measurable function h :

E k → R E[h(Z k )|Z k-1 ] = E k h(z k )Q k (dz k |Z k-1 ).
For any bounded measurable function h :

E 0 × • • • × E n → R we have: E[h(Z 0 , . . . , Z n )] = Enו••×E0 h(z 0 , . . . , z n )Q n (dz n |z n-1 ) • • • Q 1 (dz 1 |z 0 )ν 0 (dz 0 )
Let Z k = (Z 0 , Z 1 , . . . , Z k ) be a trajectory of size k, and let

E k = E 0 × E 1 × • • • × E k be
the set of trajectories of size k that we equip with the product σ-algebra

E k = E 0 ⊗ E 1 ⊗ • • • ⊗ E k .
For i < j, consider two trajectories z i and z j : when it is necessary to differentiate the coordinates of these trajectories we write the coordinates z i,k for k ≤ i and z j,k for k ≤ j such that z i = (z i,0 , z i,1 , . . . , z i,i ) and z j = (z j,0 , z j,1 , . . . , z j,j ). We introduce the Markov Chain of the trajectories (Z k ) k≥0 with values in the measurable spaces (E k , E k ), and with the transition kernels M k such that:

M k (dz k |z k-1 ) = δ z k-1 d(z k,1 , . . . , z k,k-1 ) Q k (dz k |z k-1,k-1 ).
For any bounded measurable function h : E n → R we have

p h = E [h(Z n )] = Enו••×E0 h(z n ) n k=1 M k (dz k |z k-1 )ν 0 (dz 0 ). (1)
The IPS method provides an estimator ph of p h with a different variance than the Monte-Carlo estimator. It was first introduced in [START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF], and with an alternative formulation in [START_REF] Del Moral | Sequential monte-carlo samplers[END_REF]. The IPS takes in input potential functions, denoted G k :

E k → R +
for k ≤ n. When these potential functions are properly tuned, the method can yield an estimator with a smaller variance than the Monte-Carlo estimator. When such is the case, the IPS method can assess the quantity p h with significantly less simulation runs than the Monte-Carlo method, for the same accuracy. Such variance reduction method is recommended in rare event analysis, where the Monte-Carlo method is well-known to be computationally intensive.

Rare event analyses are often carried out for reliability assessment where we want to assess the probability of failure of a reliable system. Typically there is a region D of E n that corresponds to the system failure, and we want to assess the probability of failure which happens when Z n enters D. In order to assess this probability of failure

P(Z n ∈ D), one can take h(Z n ) = 1 D (Z n ) so that p h = P(Z n ∈ D)
and use the IPS method to get the estimation ph . Although the main application of the IPS method is reliability assessment where it relates to the case h = 1 D , the result of this paper will be presented in its most general form, where h is an arbitrary measurable function.

As we said, the choice of the potential functions (G k ) k<n is paramount because it determines the variance of the IPS estimator, but so far, little information has been provided on the form of efficient potential functions. The standard approach is to find the best potential function within a set of parametric potential functions, and so the efficiency of the method strongly depends on the quality of the chosen parametric family. For instance, in [START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF] the authors obtain their best variance reduction by choosing

G k (Z k ) = exp [-λV (Z k )] exp [-λV (Z k-1 )]
where λ is a positive tuning parameter, and the quantity V (z) = a-z roughly measures the proximity of z to the critical region that was D = [a; +∞). In [START_REF] Wouters | Rare event computation in deterministic chaotic systems using genealogical particle analysis[END_REF] they stress out that it seems better to take a time-dependent proximity function V k instead of V , yielding:

G k (Z k ) = exp [-λV k (Z k )] exp [-λV k-1 (Z k-1 )]
,

where the quantities V k (z) are again measuring the proximity of z to D. Once the set of parametric potential functions is chosen, it is necessary to optimize the tuning parameters of the potentials. Different methods have been proposed. In [START_REF] Jacquemart | Tuning of adaptive interacting particle system for rare event probability estimation[END_REF], an empirical heuristic algorithm is provided; in [START_REF] Balesdent | Optimisation of interacting particle systems for rare event estimation[END_REF] a meta model of the variance is minimized; in [START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF] the large deviation principle is used as a guide. One other common option for the potential functions is the one done in splitting methods. Indeed the splitting method can also be seen as a version of the IPS method [START_REF] Cérou | Limit theorems for the multilevel splitting algorithm in the simulation of rare events[END_REF]. In this method one wants to assess the probability that a random variable Z belongs to a subset B n . A succession of nested sets

E = B 0 ⊇ B 1 ⊇ B 2 ⊇ • • • ⊇ B n is chosen by the practitioner
or possibly chosen in an adaptive manner [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF]. One considers a sequence ofrandom variables (Z i ) i=1..n such that the small probability P(Z ∈ B n ) can be decomposed into a product of conditional probabilities:

P(Z ∈ B n ) = n i=1 P(Z i ∈ B i |Z i-1 ∈ B i-1 )
and

P(Z ∈ B n ) = E [h(Z n )] by setting h(Z n ) = 1 Bn (Z n ).
In this method the potential functions are chosen of the following form

G k (Z k ) = 1 B k (Z k ).
One usually optimizes the variance reduction within this family of potential functions by optimizing the choice of the sets (B k ) k≤n .

In this paper we tackle the issue of the choice of the potential functions. Our contribution is to provide the expressions of the theoretical optimal potential functions that minimize the variance of the estimator of the IPS method. We hope these expressions will lead the practitioners to design more efficient potential functions, that are closer from the optimal ones.

The rest of the paper is organized as follows. The section 2 introduces the IPS method, and section 3 presents the potential functions that minimize the asymptotic variance of the IPS estimator, then the section 4 presents two example of applications on a toy model, and finally section 5 discusses the implications of our results.

In the rest of the paper we use the following notations: We denote by M(A) the set of bounded measurable functions on a set A. If f is a bounded measurable function, and η is a measure we note η(f ) = f dη. If M is a Markovian kernel, we denote by M (f ) the function such that M (f )(x) = f (y)M (dy|x), and for a measure η, we denote by ηM the measure such that

ηM (f ) = f (y)M (dy|x)η(dx).
2. The IPS method

A Feynman-Kac model

The IPS method relies on a Feynman-Kac model [START_REF] Moral | Genealogical and Interactig Particle Systems With Applications[END_REF] which is defined in this subsection. For each k ≤ n, we define a target probability measure ηk on (E k , E k ), such that:

∀B ∈ E k , ηk (B) = E 1 B (Z k ) k s=0 G s (Z s ) E k s=0 G s (Z s ) . ( 2
)
For each k, 1 ≤ k ≤ n, we define the propagated target probability measure η k on (E k , E k ) such that η k = ηk-1 M k-1 and η 0 = η0 . We have:

∀B ∈ E k+1 , η k+1 (B) = E 1 B (Z k+1 ) k s=0 G s (Z s ) E k s=0 G s (Z s ) . (3)
Let Ψ k be the application that transforms a measure η defined on E k into a measure Ψ k (η) defined on E k and such that

Ψ k (η)(f ) = G k (z)f (z)η(dz) η(G k ) . (4)
We say Ψ k (η) gives the selection of η through the potential G k . Notice that ηk is the selection of η k as ηk = Ψ k (η k ). The target distributions can therefore be built according to the following pattern of successive selection and propagation steps:

η k Ψ k -----→ ηk .M k -----→ η k+1 .
We also define the associated unnormalized measures γk and γ k+1 , such that for

f ∈ M(E k ): γk (f ) = E f (Z k ) k s=0 G s (Z s ) and ηk (f ) = γk (f ) γk (1) , (5)
and for f ∈ M(E k+1 ):

γ k+1 (f ) = E f (Z k+1 ) k s=0 G s (Z s ) and η k+1 (f ) = γ k+1 (f ) γ k+1 (1) . (6)
Denoting

f h (z) = h(z) n-1
s=0 Gs(z) , notice that we have:

p h = γ n (f h ) = η n (f h ) n-1 k=0 η k G k . (7)

The IPS's algorithm and its estimator

The IPS method provides an algorithm to generate samples whose weighted empirical measures approximate the probability measures η k and ηk respectively for each step k. These approximations are then used to provide an estimator of p h . For the sample approximating η k , we denote Z j k the j th trajectory and W j k its weight. Respectively in the sample approximating ηk , we denote ∼ Z j k the j th trajectory and ∼ W j k its associated weight. For simplicity reasons, in this paper, we consider that the samples all contain N trajectories, but it is possible to modify the sample size at each step, as illustrated in [START_REF] Lee | Variance estimation in the particle filter[END_REF]. The empirical measure approximating η k and ηk are denoted by η N k and ηN k and are defined by:

ηN k = N i=1 ∼ W i k δ ∼ Z i k and η N k = N i=1 W i k δ Z i k . (8)
So for all k ≤ n and

f ∈ M(E k ), ηN k (f ) = N i=1 ∼ W i k f ∼ Z i k and η N k (f ) = N i=1 W i k f Z i k . (9)
By plugging these estimations into equations ( 5) and ( 6), we get estimations for the unnormalized distributions. Denoting by γN k and γ N k these estimations, for all k ≤ n and f ∈ M(E k ), we have:

γN k (f ) = ηN k (f ) k-1 s=0 η N s (G s ) and γ N k (f ) = η N k (f ) k-1 s=0 η N s (G s ). (10)
In particular if we apply [START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF] to the test function f h , we get an estimator ph of p h defined by:

ph = η N n (f h ) n-1 k=0 η N k G k . ( 11
)
The IPS algorithm builds the samples sequentially, alternating between a selection step and a propagation step.

The k th selection step transforms the sample (Z j k , W j k ) j≤N , into the sample (

∼ Z j k , ∼ W j k ) j≤N .
This transformation is done with a multinomial resampling scheme. This means that the ∼ Z j k 's are drawn with replacement from the sample (Z j k ) j≤N , each

Initialization: k = 0, ∀j = 1..N, Z j 0 i.i.d.
∼ η 0 and W j 0 = 1 N , and

∼ W j 0 = G0(Z j 0 ) s G0(Z s 0 )
while k < n do Selection:

( Ñ j k ) j=1..N ∼ M ult N, ( ∼ W j k ) j=1..N ∀j := 1..N, ∼ W j k := 1 N Propagation:
for j := 1..N do using the kernel M k , continue the trajectory

∼ Z j k to get Z j k+1 set W j k+1 = ∼ W j k and ∼ W j k+1 = W j k+1 G k+1 (Z j k+1 ) s W s k+1 G k+1 (Z s k+1 ) if ∀j, ∼ W j k+1 = 0 then ∀q > k, set η N q = ηN q = 0 and Stop else k := k + 1
W j k G k (Z j k ) N i=1 W i k G k (Z i k )
to be drawn each time. We let Ñ j k be the number of times the particle Z j k is replicated in the sample (

∼ Z j k , ∼ W j k ) j , so N = N j=1 Ñ j k .
After this resampling the weights ∼ W j k are set to 1 N . The interest of this selection by resampling is that it discards low potential trajectories and replicates high potential trajectories. Thus, the selected sample focuses on trajectories that will have a greater impact on the estimations of the next distributions once extended.

If one specifies potential functions that are not positive everywhere, there can be a possibility that at a step k we get ∀j, G k (Z j k ) = 0. When this is the case, the probability for resampling can not be defined, the algorithm stops, and we consider that ∀s ≥ k the measures ηN s and η N s+1 are equal to the null measure.

Then the k th propagation step transforms the sample (

∼ Z j k , ∼ W j k ) j≤N , into the sample (Z j k+1 , W j k+1 ) j≤N .
Each trajectory Z j k+1 is obtained by extending the trajectory ∼ Z j k on step further using the transition kernel M k . The weights satisfy W j k+1 = ∼ W j k , ∀j. Then the procedure is iterated until the step n. The full algorithm to build the samples is displayed in table 1.

For k < n, we denote by Êk = {z k ∈ E k , G k (z k ) > 0} the support of G k , and we denote Ên = {z n ∈ E n , h(z n ) > 0} the support of h. We will make the following assumption on the potential functions:

∃ ε > 0, ∀k ≤ n, ∀z k-1 ∈ Êk-1 , M k-1 ( Êk |z k-1 ) > ε, (G)
Theorem 1. When the potential functions satisfy (G), ph is unbiased and strongly consistent.

The proof of theorem 1 can be found in [START_REF] Moral | Genealogical and Interactig Particle Systems With Applications[END_REF] chapter 7.

Theorem 2. When the potential are strictly positive:

√ N (p h -p h ) d -→ N →∞ N 0, σ 2 IP S,G (12)
where, with the convention that

-1 i=0 G i (Z i ) = -1 i=0 G -1 i (Z i ) = 1: σ 2 IP S,G = n k=0 E k-1 i=0 G i (Z i ) E E[h(Z n )|Z k ] 2 k-1 s=0 G -1 s (Z s ) -p 2 h . (13)
A proof of this CLT can be found in [START_REF] Moral | Genealogical and Interactig Particle Systems With Applications[END_REF] chapter 9. This CLT is an important result of the particle filter literature. The non-asymptotic fluctuations of the particle filters as been studied in [START_REF] Cérou | A nonasymptotic theorem for unnormalized feynman-kac particle models[END_REF]. Recently two weakly consistent estimators of the asymptotic variance σ 2 IP S,G , that are based on a single run of the method, have been proposed in [START_REF] Lee | Variance estimation in the particle filter[END_REF]. One of these estimators is closely related to what was proposed in [START_REF] Peng Chan | A general theory of particle filters in hidden markov models and some applications[END_REF].

The theoretical expression for the optimal potential

Here we aim at estimating p h , the finality is not the estimation of some distributions η k and ∼ η k for k < n, so we can choose any potential functions (G k ) k<n providing they are positive. But the choice of the potential functions has an impact on the variance of the estimation, so we would like to find potential functions that minimize the asymptotic variance [START_REF] Jacquemart | Tuning of adaptive interacting particle system for rare event probability estimation[END_REF]. Also, note that if potential functions G k and G k are such that G k = a.G k with a > 0, then they yield the same variance: σ 2 IP S,G = σ 2 IP S,G . Therefore all potential functions will be defined up to a multiplicative term. Theorem 3. For k ≥ 1, let G * k be defined by:

G * k (Z k ) ∝            E E h(Zn) Z k+1 2 Z k E E h(Zn) Z k 2 Z k-1 if E E h(Z n ) Z k 2 Z k-1 = 0 0 if E E h(Z n ) Z k 2 Z k-1 = 0 ( 14
)
and for k = 0,

G * 0 (Z 0 ) ∝ E E h(Z n ) Z 1 2 Z 0 . ( 15
)
The potential functions minimizing σ 2 IP S,G are the ones that are proportional to the G * k 's ∀k ≤ n. The optimal variance of the IPS method with n steps is then

σ 2 IP S,G * = E E h(Z n ) Z 0 2 -p 2 h + n k=1    E E E h(Z n ) Z k 2 Z k-1 2 -p 2 h    . (16)
Proof. As we lack mathematical tools to minimize σ 2 IP S,G over the set of positive functions (G k ) k≤n , we had to guess the expressions ( 14) and ( 15) before providing the proof of the results. We begin this proof by presenting the heuristic reasoning that provided the expressions ( 14) and [START_REF] Wouters | Rare event computation in deterministic chaotic systems using genealogical particle analysis[END_REF].

Assuming we already know the k -2 first potential functions, we started by trying to find the k -1-th potential function G k-1 that minimizes the k-th term of the sum in [START_REF] Jacquemart | Tuning of adaptive interacting particle system for rare event probability estimation[END_REF]. This is equivalent to minimize the quantity

E k-1 i=0 G i (Z i ) E E E[h(Z n )|Z k ] 2 Z k-1 k-1 s=0 G -1 s (Z s ) (17) over G k-1 .
As the G k-1 are equivalent up to a multiplicative constant, we simplify the equation by choosing a multiplicative constant so that

E k-1 i=0 G i (Z i ) = 1.
Our minimizing problem then becomes the minimization of (17) under the constraint

E k-1 i=0 G i (Z i ) = 1.
In order to be able to use a Lagrangian minimization we temporarily assume that the distribution of Z k-1 is discrete and that Z k-1 takes its values in a finite or numerble set E. For z ∈ E, we denote a z = P(Z k-1 = z) and

d z = E E[h(Z n )|Z k ] 2 Z k-1 = z and g z = k-2 i=0 G i (Z i )G k-1 (z) our minimization problem becomes the minimization of L = z∈E p z d z g z -λ 1 - z∈E p z g z (18)
Finding the minimum of this Lagrangian we get that

g z = √ dz z ∈E p z √ d z
. Now relaxing the constraint of the multiplicative constant, we get that

k-1 i=0 G i (Z i ) ∝ E E[h(Z n )|Z k ] 2 Z k-1 ,
which gives the desired expressions. After these heuristic arguments we can now rigorously check that these expressions, obtained by minimizing each of the term of sum in [START_REF] Jacquemart | Tuning of adaptive interacting particle system for rare event probability estimation[END_REF] one by one, also minimize the whole sum for any distribution of the Z k-1 's.

The proof now consists in showing that, for any set of potential functions (G s ), we have σ 2 IP S,G ≥ σ 2 IP S,G * . This is done by bounding from below each term of the sum in [START_REF] Jacquemart | Tuning of adaptive interacting particle system for rare event probability estimation[END_REF].

We start by decomposing a product of potential functions as follows:

∀k ∈ {1, . . . , n}, k-1 s=0 G s (Z s ) = k-1 (Z k-1 ) k-1 s=0 G * s (Z s ) + ¯ k-1 (Z k-1 ) (19)
where when

Z k-1 ∈ supp k-1 s=0 G * s , k-1 (Z k-1 ) = k-1 s=0 G s (Z s ) k-1 s=0 G * s (Z s )
, and ¯ k-1 (Z k-1 ) = 0 and when

Z k-1 / ∈ supp k-1 s=0 G * s , k-1 (Z k-1 ) = 0, and ¯ k-1 (Z k-1 ) = k-1 s=0 G s (Z s ).
Using (19) we get that

E k-1 s=0 G s (Z s ) E E[h(Z n )|Z k ] 2 k-1 s=0 G -1 s (Z s ) = E k-1 (Z k-1 ) k-1 s=0 G * s (Z s ) E E E[h(Z n )|Z k ] 2 Z k-1 k-1 s=0 G -1 s (Z s ) + E ¯ k-1 (Z k-1 ) E E E[h(Z n )|Z k ] 2 Z k-1 k-1 s=0 G -1 s (Z s ) ≥ E k-1 (Z k-1 ) k-1 s=0 G * s (Z s ) E E E[h(Z n )|Z k ] 2 Z k-1 k-1 s=0 G s (Z s ) + 0 (20) For Z k-1 ∈ supp k-1 s=0 G * s we have: k-1 s=0 G * s (Z s ) ∝ E E h(Z n ) Z k 2 Z k-1 So supp E E[h(Z n )|Z k ] 2 Z k-1 = supp k-1 s=0 G * s and we get E E E[h(Z n )|Z k ] 2 Z k-1 k-1 s=0 G s (Z s ) = E E E[h(Z n )|Z k ] 2 Z k-1 k-1 (Z k-1 ) k-1 s=0 G * s (Z s ) = E 1 k-1 (Z k-1 ) k-1 s=0 G * s (Z s ) . (21)
Combining ( 21) with inequality (20) we get that

E k-1 s=0 G s (Z s ) E E[h(Z n )|Z k ] 2 k-1 s=0 G -1 s (Z s ) ≥ E k-1 (Z k-1 ) k-1 s=0 G * s (Z s ) E 1 k-1 (Z k-1 ) k-1 s=0 G * s (Z s ) (22)
and using the Cauchy-Schwarz inequality on the right term, we get that

E k-1 s=0 G s (Z s ) E E[h(Z n )|Z k ] 2 k-1 s=0 G -1 s (Z s ) ≥ E k-1 s=0 G * s (Z s ) 2 = E k-1 s=0 G * s (Z s ) E E E[h(Z n )|Z k ] 2 Z k-1 k-1 s=0 G * s (Z s ) . ( 23
)
By summing the inequalities (23) for each k, we easily see that

σ 2 IP S,G ≥ σ 2 IP S,G * ,
which completes the proof of the theorem.

Remark that the optimal potential is unfortunately not positive everywhere, so it violates the hypothesis under which the TCL was proven in [START_REF] Moral | Genealogical and Interactig Particle Systems With Applications[END_REF]. We claim that this question of the positiveness of the potential has not much interest in practice. Assume we take potential functions (G k ) k<n such that in the equation (20) we have k (Z k ) = 1 and ¯ k (Z k ) = ε > 0 with ε very small. Choosing ε small enough, we can get (G k ) k<n as close has we want from (G * k ) k<n . With such potential functions (G k ) k<n and ε small enough, it is very likely that we would obtained the same samples in the algorithm as if we had taken the potentials (G * k ) k<n , and so we would have the same estimation.

Moreover, with such potential functions we would also have a TCL with a variance very close to σ 2 IP S,G * . (According to (20) by choosing ε as small as we want, we get a variance as close as we want from σ 2 IP S,G * .) In practice, positive potential functions (G k ) k<n with ε close to zero give the same results as the potentials (G * k ) k<n .

Application on a toy model

In this section we apply the IPS method on a toy system for which we have explicit formulas. The system under consideration is the Gaussian random walk Z k+1 = Z k + ε k+1 , Z 0 = 0, where the (ε k ) k∈{1,...,n} are i.i.d. Gaussian random variables with mean zero and variance one. We explore two situations, one where the quantity to estimate,the optimal potential and the variance of the estimator can be calculated explicitly, and one where these quantities can be approximated by a large deviation inequality.

First example

In the first situation, taking b, a > 0 and n ∈ N\{0}, the goal is to compute the

expectation p h = E[h(Z n )] when h(z) = exp[b(z -a)].
As Z n is a centered Gaussian of variance n, a simple calculation gives that for k < n:

E[h(Z n )|Z k = z] = exp (n -k) b 2 2 + b(z -a) . (24)
Consequently we have that:

p h = exp n 2 b 2 -ab , (25)
and that, for k ≥ 1:

E E h(Z n ) Z k+1 2 Z k E E h(Z n ) Z k 2 Z k-1 = exp - b 2 2 + b(Z k -Z k-1 ) (26)
with

E E h(Z n ) Z 1 2 Z 0 = exp (n + 1)b 2 2 + b(Z 0 -a) , (27)
or equivalently that:

G * k (Z k ) ∝ exp b(Z k -Z k-1 ) (28)
with

G * 0 (Z 0 ) ∝ exp bZ 0 . (29)
Using the equations ( 28) and ( 29) in [START_REF] Jacquemart | Tuning of adaptive interacting particle system for rare event probability estimation[END_REF], it can easily be shown that the variance of the IPS estimator with these optimal potential functions is:

σ 2 IP S,G * = n exp(b 2 ) -1 exp nb 2 -2ab = n exp(b 2 ) -1 p 2 h . (30)
It is notable that we obtain an asymptotically optimal variance [START_REF] Dupuis | Importance sampling, large deviations, and differential games[END_REF], i.e. a variance proportional to p 2 h . In order to confirm these theoretical results, we have carried out a simulation study.

We have run the method 200 times with N = 10 5 , n = 10 and different values of a and b, and for each of these values we have computed the mean of the estimation and the empirical variance of the estimation. The results are displayed in table 2, where we compare the theoretical value of p h to the empirical mean of the 200 estimations, and

σ 2
IP S,G * to the empirical variance. As the theoretical values are close to the empirical ones, this confirms that the method is unbiased and that the variance given by the equation ( 13) is the right one. We also compare the empirical variance to the variance of the Monte-Carlo estimator, showing that, on this example, the IPS method provides a significant variance reduction with the optimal potential, as the variance is reduced by at least a factor 10 4 on the considered cases.

Second example

In the second situation, the goal is to compute the probability that Z n exceeds a large positive value a. Therefore we take h(

Z n ) = 1 [a;+∞) (Z n ) so that p h = P(Z n ≥ a). (31)
In that case one can not compute E[h(Z n )|Z k = z] but the Chernov-Bernstein's inequality gives the following sharp exponential bound:

E[h(Z n )|Z k = z] ≤ exp - (a -z) 2 2(n -k) , (32)
from which we can deduce that:

k-1 i=0 G * i (Z i ) = E E h(Z n ) Z k 2 Z k-1 ≤ C 1 exp C 2 - (Z k-1 -a) 2 2(n -k + 2) , (33)
N = 2 * 10 5 , n = 10
where C 1 and C 2 are some constants independent of Z k-1 . One can therefore try to set the product of the potentials equal to this upper bound, which yields:

for k ≥ 1, G k (Z k ) ∝ exp - (Z k -a) 2 2(n -k + 1) + (Z k-1 -a) 2 2(n -k + 2) (34) and G 0 (Z 0 ) ∝ exp - (Z 0 -a) 2 2(n -1) . ( 35
)
Similarly as for the first example, we have carried out a simulation study. We have run the method 200 times with the potentials defined in equation ( 34) and ( 35)

and taking N = 2 × 10 5 , and n = 10, and different values of a, and for each of these values we have computed the empirical mean of the estimation and the empirical variance of the estimation. The results are displayed in table 3. We compare these estimations with the actual values of p h and the variance of the Monte-Carlo method

σ 2 M C
, showing that the potentials built with the Chernov-Bernstein large deviation inequality and our formula yield a significant variance reduction. Indeed the variance reduction compared to the Monte-Carlo method is at least by a factor 8500, and at best by a factor 2.4 * 10 6 . We also compared the efficiency of different potentials. results obtained with N = 10 5 and n = 10

a p h σ 2 IP S,G * σ 2 M C mean(p h ) σ2 IP S,G 4
We run the method with 1) the potential used on a Gaussian random walk in [START_REF] Jacquemart | Tuning of adaptive interacting particle system for rare event probability estimation[END_REF]:

G k (Z k ) = exp [α(Z k -Z k-1 )
] where the parameter was optimized to α = 1.1, 2) the potential built with the Chernov-Bernstein large deviation inequality, and 3) with the optimal potential that we computed using Gaussian quadrature formulas. The results are displayed in table 4, and show that indeed the potential functions (G * k ) k<n , where results obtained for p h = 1.05 * 10 -6 N = 2000, n = 10, a = 15, α = 1.1

G * k (Z k ) = R ∞ a exp - (zn-z k+1 ) 2 2(n-k+1) dzn 2 exp - (z k+1 -Z k ) 2 2 dz k+1 R ∞ a exp - (zn -z k ) 2 2(n-k) dzn 2 exp - (z k -Z k-1 ) 2 2 dz k , yield the best variance. G k (Z k ) mean(p h ) σ2 IP S,G exp [α(Z k -Z k-1 )] 1.04 * 10 -6 2.44 * 10 -10 exp -(Z k -a) 2 2(n-k+1) + (Z k-1 -a) 2 2(n-k+2)

Discussion and implication

In this paper we give a closed form expression of the optimal potential functions for the IPS method with multinomial resampling, and for its minimal variance. The existence of optimal potential function, proves that the possible variance reduction of an IPS method is lower-bounded. The expressions have been validated analytically, and the expression for the minimal variance has been empirically confirmed in toy examples.

Furthermore the results found in the literature seem consistent with our findings. Indeed, in [START_REF] Del Moral | Genealogical particle analysis of rare events[END_REF] the authors made the observation that it seemed better to build a potential which increments of an energy function, this observation is confirmed as the optimal potential is the multiplicative increment of the quantity

E E h(Z n ) Z k 2 Z k-1
which is then the optimal energy function. Also, the fact that in [START_REF] Wouters | Rare event computation in deterministic chaotic systems using genealogical particle analysis[END_REF] the authors find better results with time-dependent potential is explained by the fact the expression of the optimal potential shows a dependency on k. Finally, as splitting methods can be viewed as a version of the IPS-method with indicator potential functions, our results

show that the selections of splitting algorithms are not optimal, and could be improved by using information on the expectations E h(Z n ) Z k = z .

The optimal potential functions may be hard to find in practice. Indeed, the expectations E h(Z n ) Z k = z play a big role in the expression of the optimal potentials, but if we are trying to assess p h = E h(Z n ) , we typically lack information about the expectations E h(Z n ) Z k = z . If no information on the expectation E h(Z n ) Z k = z is available, it might be preferable to use more naive variance reduction method, where no input functions are needed. In such context, the Weighted Ensemble (WE) method [START_REF] Aristoff | Analysis and optimization of weighted ensemble sampling[END_REF][START_REF] Aristoff | Optimizing weighted ensemble sampling for steady states and rare events[END_REF] seems to be a good candidate, as it does not take in input potential functions but only a partition of the state space. Conversely if the practitioner has information about the expectation E h(Z n ) Z k = z , this information could be used to derive very efficient potentials.

The knowledge of these expectations is therefore crucial for a well optimized use of the IPS method, but it is interesting to remark that the same knowledge seem to be crucial for a well optimized importance sampling method. Indeed the optimal density of importance sampling (that gives a zero variance) depends on

E h(Z n ) Z [0,t] = z [0,t]
when it is used on a piecewise deterministic process [START_REF] Chraibi | Speed-up reliability assessment for multi-component systems: importance sampling adapted to piecewise deterministic markovian processes[END_REF]. This confirms the well known fact that, with a good knowledge of the dynamic of the process (Z k ) k≥0 , the importance sampling methods is preferable to the IPS. Nonetheless the IPS method may still be preferred to the importance sampling methods when that knowledge is not available.

When one fears to be in an over-biasing situation the IPS may be preferable, as in the IPS method we do not alter the propagation, the over-biasing phenomenon should then be less important than with importance sampling methods.

Table 1 :

 1 IPS algorithm

	trajectory Z j k having a probability

  ) 4.50 × 10 -10 8.13 × 10 -19 4.50 × 10 -10 4.50 × 10 -10 7.39 × 10 -17

	35	35	35	35	40	40	40	40	a
	log(1 + 4 n ) 8.19 × 10 -9	log(1 + 3 n ) 6.07 × 10 -8	log(1 + 2 n ) 8.04 × 10 -7	log(1 + 1 n ) 3.27 × 10 -5	log(1 + 4	log(1 + 3 n ) 4.69 × 10 -9	log(1 + 2 n ) 9.51 × 10 -8	log(1 + 1 n ) 6.98 × 10 -6	b p h
	2.69 × 10 -16 8.19 × 10 -9 8.19 × 10 -9	1.11 × 10 -14 6.07 × 10 -8 6.07 × 10 -8	1.29 × 10 -12 8.04 × 10 -7 8.04 × 10 -7	1.09 × 10 -9 3.27 × 10 -5 3.27 × 10 -5		6.61 × 10 -17 4.69 × 10 -9 4.69 × 10 -9	1.81 × 10 -14 9.51 × 10 -8 9.51 × 10 -8	6.87 × 10 -16 6.98 × 10 -6 6.98 × 10 -6	σ 2 IP S,G * σ 2 M C mean(p h )
	2.85 × 10 -16	1.06 × 10 -14	1.35 × 10 -12	9.69 × 10 -10		7.37 × 10 -17	1.76 × 10 -14	6.72 × 10 -16	σ2 IP S,G *

n

Table 2 :

 2 Theoretical and empirical comparisons (example 1)

Table 3 :

 3 Theoretical and empirical comparisons (example 2)

Table 4 :

 4 Comparisons of the efficiency of potentials (example 2)