Timothé Collet Centralesupélec

Olivier Pietquin
email: olivier.pietquin@univ-lille1.fr

Bayesian Credible Intervals for Online and Active Learning of Classification Trees

Classification trees have been extensively studied for decades. In the online learning scenario, a whole class of algorithms for decision trees has been introduced, called incremental decision trees. In the case where subtrees may not be discarded, an incremental decision tree can be seen as a sequential decision process, consisting in deciding to extend the existing tree or not. This problem involves an trade-off between exploration and exploitation, which is addressed in recent work with the use of Hoeffding's bounds. This paper proposes to use Bayesian Credible Intervals instead, in order to get the most out of the knowledge of the output's distribution's shape. It also studies the case of Active Learning in such a tree following the Optimism in the Face of Uncertainty paradigm. Two novel algorithms are introduced for the online and active learning problems. Evaluations on realworld datasets show that these algorithms compare positively to state-of-the-art.

I. INTRODUCTION

We place this paper in the context of classification. This supervised learning framework consists in learning a mapping from an instance space to a finite discrete label set. For this purpose, it uses a training set containing a finite number of instances labeled by an expert. Two scenarios exist concerning the acquisition of data. In the batch learning scenario, the entire training set is given at once, and the classifier is computed considering all the instance/label pairs at the same time. In the online learning scenario, the data become available only sequentially, with the current classifier being required at every time step. This is the case, for example, in a system which learns while being used, with direct feedback. In this second scenario, a classifier capable of incremental updates is preferred to one being recomputed from scratch every time it receives new data, mainly for its lower computational time.

Our work focuses on classification trees. They are acyclic directed graphs whose nodes contain a test on instances, and branches represent the outcome of this test. Each leaf is associated to a label prediction based on the distribution of labels it received. In order to maximize the performance, the tree aims to partition the instance space into homogeneous labels zones. Learning this classifier consists in reducing a given measure of heterogeneity, e.g. information entropy or standard deviation, estimated from received labels in the leaves, by performing recursive splits. An example is the c4.5 algorithm [START_REF] Quinlan | C4. 5: programs for machine learning[END_REF], it uses the information gain, or difference of entropy, to define the best splits.

Classification trees in the online learning scenario are called online or incremental classification trees [START_REF] Crawford | Extensions to the cart algorithm[END_REF] [3] [START_REF] Domingos | Mining high-speed data streams[END_REF].

We are only interested in the case where subtrees cannot be discarded and the existent tree can only be extended. This way, an online classification tree dynamically builds the tree, which corresponds to adaptively partition the instance space. It is a sequential decision process [START_REF] Degroot | Optimal statistical decisions[END_REF] which, for each newly received instance-label pair, considers changing the concerned leaf into one or more decision nodes. Once a test has been installed, it cannot be removed, it is thus critical to be sure that it is needed before installing it.

Learning an online classification tree involves an exploration/exploitation dilemma. Indeed, ideally, a split is needed if it decreases the true heterogeneities of the leaves instead of their estimates based on received labels. When considering to install a decision node leading to a decrease of the estimates, we can either decide to effectively split, which results in exploitation of the estimates. Or we can decide to wait for a better estimate, which results in exploration.

The exploration/exploitation dilemma has been successfully addressed in the case of search trees in [START_REF] Kocsis | Bandit based monte-carlo planning[END_REF], which introduces the UCT algorithm. In this problem, each leaf is associated to a distribution of rewards and the algorithm aims to find the best leaf through repeated playouts. Based on Upper Confidence Bounds [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF], it builds confidence intervals on expected rewards and follows the branch for which the upper bound is highest. In the case of online classification trees, only one playout is performed. Thus, only trustworthy actions must be taken by considering lower bounds instead of upper bounds. In [START_REF] Domingos | Mining high-speed data streams[END_REF] and [START_REF] Gama | Forest trees for on-line data[END_REF], the authors use Hoeffding's bounds on the gain of heterogeneities, and split if the lower bound is positive. Their experiments show good results on relatively large datasets. However, we show that for more limited datasets, those algorithms are not very efficient. Indeed, Hoeffding's bounds are known to be inaccurate with a small number of samples.

If the family of an output distribution is known a priori, it is possible to use the family of conjugate priors along with Bayesian inference to build confidence intervals. These intervals are known as Bayesian Credible Intervals. In [START_REF] Kaufmann | On bayesian upper con-fidence bounds for bandit problems[END_REF], the authors show that using such bounds prove optimal for the multi-armed bandit. In addition, experiments show that their performances are better than when using Hoeffding's bounds.

In this paper, we introduce an online classification tree algorithm using Bayesian Credible Intervals. Indeed, in the context of classification, the family of the output distribution is Bernoulli. Bayesian Credible Intervals can therefore be used. As they come from concentration equalities instead of inequalities, they are the tightest confidence bounds possible whatever the number of sample, and are thus well fitted for online learning.

We also study the particular case of Active Learning [START_REF] Settles | Active learning literature survey[END_REF]. An Active Learning algorithm is used to decide which instances to label. Indeed, the use of an expert for annotation is considered time-and/or money-consuming. Given a finite budget of annotations, an Active Learning algorithm aims to maximize the performance of the classifier. Although some one-shot Active Learning algorithms exist [11] [12], in which all the instances to be labeled are chosen at once, Active Learning is usually seen as a sequential process [START_REF] Degroot | Optimal statistical decisions[END_REF]. Indeed, the choice of each instance is based on the current predictions, depending on all the labels received so far. An online classifier is thus necessary.

At every moment, the online classification tree defines a partition of the instance space. An Active Learning strategy thus selects the leaf from which to draw the instance. In [START_REF] Collet | Active learning for classification: An optimistic approach[END_REF], the authors introduce an Active Learning algorithm for the case of single partitions. Although, the partition used cannot change all along the process of the algorithm. Thus, the partition cannot adapt to the distribution of labels among the instance space. Moreover, the number of subsets in the partition does not depend on the number of labeled instances. This makes it impractical to use on real world datasets.

We propose to adapt their selection strategy to online classification trees. This suppose that the instances were taken at random inside each leaf. By not allowing to remove any installed tests, the online classification tree respects this constraint.

This paper is structured as follows. In section III, we introduce the Bayesian Online Classification Tree (BOCT) algorithm. We thus show how to compute Bayesian Credible Intervals and use them in place of Hoeffding bounds. In section IV, we introduce an extension to BOCT including Active Learning. We derive our own selection strategy inspired by [START_REF] Collet | Active learning for classification: An optimistic approach[END_REF]. In section V, we present the experiments used to evaluate the different algorithms. In section VI, we describe the results of the evaluation and compare to state-of-the-art algorithms. We show that there is a significant improvement brought by using Bayesian Credible Intervals instead of Hoeffding bounds, and that Active Learning could further improve the performance. We conclude in section VII.

II. RELATED WORK

In [START_REF] Carpentier | De l'échantillonage optimal en grande et petite dimension[END_REF], the author introduce an algorithm called Deep-MC-UCB. They study the problem of Active Learning in an adaptive partition, however with some differences. First, the purpose is to estimate uniformly well the mean values of the output distribution in each subset. Thus, it corresponds to regression rather than classification. In [START_REF] Collet | Active learning for classification: An optimistic approach[END_REF], the authors showed that using a loss specifically adapted for classification showed a significant improvement. Second, they use a predefined maximum partition, which is hard to use in a realworld problem. Yet, this motivated our toy problem. Last, they use Bernstein's bounds, where we can use Bayesian Credible Intervals, due to the fact that we study classification.

III. ONLINE CLASSIFICATION TREE

In this section, we introduce an algorithm for learning a classification tree in an online fashion. Here, we only consider binary trees. In such a tree, each non-leaf node is associated to a test on an attribute, and branches represent the outcome of the test. Thus, the nodes correspond to subsets of the instance space conditioned by the successive tests. Each leaf is given a label which is the prediction for an instance belonging to the corresponding subset.

Learning this tree corresponds to determine which test is best to install, on which attribute, and whether to install it or not. And also which label to give to each leaf. In the online learning scenario, updating the tree corresponds to adding more decision nodes to the existing tree, without the possibility to remove any. Thus, a test must be carefully examined before being permanently installed. The solution is to first study in which case would the split be useful. Then, to study the probability of such a case, conditioned to the observations. The algorithm is reported in Fig. 2 and Fig. 3. Apart from problem dependent parameters, it takes one parameter as input: δ, which is the probability of the true heterogeneity being inside the established intervals. The rest of the section describes the algorithm.

The labels drawn from one particular subset follows a Bernoulli distribution of unknown parameter. The true heterogeneity of a subset is a function of this parameter, usually information entropy, but other exist such as standard deviation or variance. A split is useful if it makes the heterogeneity decreases inside each obtained leaf. Thus, if the weighted sum of true heterogeneities of children nodes is lower than the true heterogeneity of the parent node. The true heterogeneity is not known but Bayesian Credible Intervals can be computed using labels received so far. Those intervals contain the true heterogeneity with a given probability. If the lower bound of the parent node is higher than the upper bound of the children's nodes, it is highly probable the test must be installed. Fig. 1 shows the Bayesian Credible Intervals on the true heterogeneities of the parent node compared to ones on the weighted sum of true heterogeneities after several tests. In the following, we start by showing how to compute a Bayesian Credible Interval for the true heterogeneity. In this paper, we use the following notations. Suppose that the instance space is X and the label set is Y . We are only interested in binary classification, in which the label set contains two elements, in our case Y = {0, 1}. At every time step t ∈ 1, n , the algorithm receives a new instance x t ∈ X labeled y t ∈ Y and update the classification tree.

Let us consider a generic subset X s ⊂ X, we denote by T s,t the count of received instances that lie in that subset:

T s,t = t i=1 1 xi∈Xs (1)
We denote by μs,t the average of labels received in subset X s :

μs,t = 1 T s,t t i=1 1 xi∈Xs y i (2)
Let us denote µ s the parameter of the Bernoulli distribution associated to subset X s . Note that the family of conjugate priors for Bernoulli distributions is the Beta distributions. Thus, by Bayesian inference:

∀ s ∈ [0, 1] P(µ s ≥ s ) = I s (T s,t μs,t + 1, T s,t (1 -μs,t ) + 1),
where I (., α, β) is the cumulative distribution function of a Beta distribution of parameters α and β.

Let H(µ s ) be the true heterogeneity of subset X s . Usual true heterogeneity functions are:

• information entropy, with H(µ s ) = -µ s log µ s -(1 -µ s ) log 1 -µ s , • variance, with H(µ s ) = µ s (1 -µ s ), • standard deviation, with H(µ s ) = µ s (1 -µ s ).
But any function could work as long as it is concave, symmetrical with respect to 1 2 , and null for µ s = 0 and µ s = 1. It is thus, maximum for µ s = 1 2 . For any s ∈ (0, H( 1 2 )), let µ low ( s ) and µ upp ( s ) be such that H(µ low ) = H(µ upp ) = s and µ low < µ upp .

For example, if the heterogeneity is measured by the standard deviation, H(µ s ) = µ s (1 -µ s ). Then, the so-

lutions of H(µ s ) = s are µ low ( s ) = 1 2 -1 4 -2 s and µ upp ( s ) = 1 2 + 1 4 -2 s .
Let us define the function

f such that, ∀ s ∈ (0, H( 1 2 )) f ( s ) = P(H(µ s ) ≥ s ) (3) = I µupp( s) (T s,t μs,t + 1, T s,t (1 -μs,t ) + 1) -I µ low ( s) (T s,t μs,t + 1, T s,t (1 -μs,t ) + 1).
The function f is strictly decreasing, it thus admits an inverse f -1 . The upper bound s (δ) of the true heterogeneity according to any given probability δ: ∀δ ∈ [0, 1]

s (δ) = f -1 (δ). (4)
Note that the lower bound is s (1 -δ). Indeed, ∀δ ∈ [0, 1]

P(H(µ s ) ≤ s (1 -δ)) = 1 -P(H(µ s ) ≥ s (1 -δ)) = δ.
Thus, the Bayesian Credible Interval is such that the true heterogeneity is in

[ s (1 -δ), s (δ)] with probability 1 -2δ.
A low probability δ is usually chosen in order to get high probability confidence intervals.

Core algorithm

Require: X, (x i ) i∈ 1,n , δ Initialize Tree as a one leaf tree. This leaf is associated to the whole instance space X. for t = 1 to n do Receive a label y t for the instance x t Let node be the node of Tree which x t belongs to SubTree = U pdate(node, (x i ) i∈ 1,n , (y i ) i∈ 1,t , δ) Replace node by SubTree end for Ensure: Tree

(x i ) i∈ 1,n ,(y i ) i∈ 1,t , δ
Compute T node , μnode and W node using Eq. (1) (2) (5) Compute node (1 -δ) using Eq. (3) ( 4)

for a = 1 to N a do (x sorted i,a ) i∈ 1,t = sort((x i,a ) i∈ 1,t ) for i = 1 to t -1 do if x sorted i,a = x sorted i+1,a then cutValue i,a = x sorted i,a +x sorted i+1,a 2
Let X c1 and X c2 be two subsets resulting from splitting X node along attribute a with cut value

cutValue i,a Compute T c1 , μc1 , W c1 , T c2 , μc2 , W c2 (1) (2) (5) Compute c1 (δ) and c2 (δ) using Eq. (3) (4) c i,a = W c1 c1 (δ) + W c2 c2 (δ) end if end for end for {i best , a best } = arg min i,a c i,a if c i best ,a best < W node node (1 -δ) then
Create two branches from node with test:

x .,a best > cutValue i best ,a best leading to leaves child 1 and child In order not to overcharge the notations, it is left to the reader's understanding that the index s represents T s,t and μs,t . Now that we know how to compute Bayesian Credible Intervals, let us define the condition on which a test is installed. Let us consider the installation of a test node which splits a current leaf's subset X p ⊂ X into two subsets X c1 ⊂ X and

2 child 1 = U pdate(child 1 , (x i ) i∈ 1,n , (y i ) i∈ 1,t , δ) child 2 = U pdate(child 2 , (x i ) i∈ 1,n , (y i ) i∈ 1,t , δ) end if Ensure: Tree
X c2 ⊂ X. Note that X c1 ∪ X c2 = X p and X c1 ∩ X c2 = ∅.
The heterogeneities of each node must be weighted by the size of the node's subset. We denote by W s the size of a generic subset X s . This size could be the volume of the subset relatively to some metric. It could be easily computed by starting with a delimited instance space and knowing the splitting criterion. However, this does not consider the distribution of instances among the instance space. In order to estimate it, we could count the instances received in the subset. In this case, W s = T s,t . However, in order to get more precise estimates that does not change over the time, we place ourselves in a semi-supervised context. In this context, all the unlabeled instances are known from the beginning. The counts of instances now become:

W s = n i=1 1 xi∈Xs . ( 5
)
Nota Bene: the difference is the sum ends at n whereas before it ended at the current time t.

In the case where the true heterogeneities of the different subsets are known, the test should be installed if:

W p H(µ p ) > W c1 H(µ c1 ) + W c2 H(µ c2 ). (6)
However, the true heterogeneities are not known. But, we have access to Bayesian Credible Intervals which gives an information about their value. We now show a second condition which guaranties a certain probability that the previous condition is respected. We first give the upper bound of the Bayesian Credible Interval for the right side of Eq. ( 6).

Combining the two following upper bounds:

P(H(µ c1 ) ≥ c1 (δ)) = δ and P(H(µ c2 ) ≥ c2 (δ)) = δ,
we obtain

P(W c1 H(µ c1 )+W c2 H(µ c2 ) ≥ W c1 c1 (δ)+W c2 c2 (δ)) ≤ 2δ.
We remind that we have the following lower bound of the Bayesian Credible Interval for the left side of Eq. ( 6):

P(H(µ p ) ≤ p (1 -δ)) = δ.
Finally, if the lower bound of the left side of Eq. ( 6) is higher than the upper bound of the right side, we can assume that there is a high probability that the condition is met. Formally,

if W p p (δ) > W c1 c1 (δ) + W c2 c2 (δ), then P(W p H(µ p ) ≤ W c1 H(µ c1 ) + W c2 H(µ c2 )) ≤ 3δ.
This criterion is used to decide if a test has to be installed or not. However, when considering the split of a leaf, it happens that several tests may be valid. In this case, the best test to install is the one leading to the best decrease in true heterogeneity. Thus, we consider that the test for which the upper bound of children subsets' heterogeneity is lowest is best. Note that for a given parent node, the choice of the test is exactly the choice of {c 1 , c 2 }.

test best = arg min test W c1 c1 (δ) + W c2 c2 (δ).
At time step t, this process of choosing a test is carried out on the leaf of the current tree in which the newly received instance lies. Indeed, it is the only leaf affected by the new labeled instance. It is repeated for each newly created leaf every time a test is successfully installed.

In this section, the Bayesian Online Tree algorithm was introduced. It allows a classification tree to be updated when faced to an online scenario, which avoids to recompute the whole tree at each time step. This decreases the computational time because it does not reconsider to split leaves other than the one currently receiving an instance. Here, we suppose that the function f -1 can be store in a table. Indeed, the inputs are discrete and T s,t is bounded by n, and the same value is accessed many times during the process of the algorithm.

One other advantage of this method is that the installed tests are never removed. This is particularly interesting in Active Learning, which is the focus of the next section.

IV. ACTIVE LEARNING OF CLASSIFICATION TREES

In this section, we show how the learning rate of the previous classifier can be improved by actively selecting the instances to be labeled. At each time step t, the current classification tree partitions the instance space into N t subsets. As it makes no distinction between the instances belonging to the same subset, the selection strategy only choses at each time step the subset k t in which to pick an instance at random. The goal is here to design a criterion based on labels received so far, that represents the need of receiving more labels in a subset

The algorithm is described in Fig. 4. The input parameter δ is the probability on which are built all the confidence intervals. The rest of the section describes the operation of the algorithm.

Let us start with an intuition, if one subset only received labels of the same class, it is likely to receive the same label again. Thus, the predicted label will not change for this subset. On the contrary, if it received a mix of labels, getting another one not only refines the prediction but may also allow the subset to be split.

The aim of Active Learning is to get the best performance of the classifier given a finite budget of labeled instance. This performance is measured by the true risk, based on the 0/1 loss it represents the misclassification rate on unseen instances. Let us define a selection strategy which aims at decreasing the true risk. Let R t be the true risk of a classification tree at time step t, it is the sum of the true risk in each subset:

R t = Nt s=1 R s,t .
The true risk is thus bounded by

R t ≤ N t max s∈ 1,Nt R s,t .
Thus, in order to decrease the true risk, we aim to decrease the maximum true risk among subsets. For this purpose, the solution is to sample the subset for which the true risk is currently maximum:

k t = arg max s∈ 1,Nt R s,t .
Let us now specify the true risk of a subset s. The labels received in subset s follow a Bernoulli distribution of parameter µ k . Suppose that the average of received labels is μs,t , thus, the label [μ s,t ] is predicted for this subset, where [.] is the round operator. Then,

R s,t = W k µ k if [μ s,t ] = 0 W k (1 -µ k ) if [μ s,t ] = 1,

Active Learning algorithm

Require: X, (x i ) i∈ 1,n , δ 1 , δ 2 Initialize Tree as a one leaf tree. This leaf is associated to the whole instance space X. for t = 1 to n do for s = 1 to N t do Compute μs,t , T s,t and W k Compute k (δ 2 ) using Eq. ( 7) end for Compute k t = arg max k k Receive a label y t for a random instance x t ∈ X kt Let node be the node of Tree which x t belongs to SubTree = U pdate(node, (x i ) i∈ 1,n , (y i ) i∈ 1,t , δ) Replace node by SubTree end for Ensure: Tree Note that, even though this optimal selection strategy is inspired from [START_REF] Collet | Active learning for classification: An optimistic approach[END_REF], it is not exactly identical. In their work, the authors chose to sample according to the maximum difference between the current true risk and the optimal true risk. Indeed, in the case of a fixed partition, it is unnecessary to spend effort to improve the prediction if the risk cannot be decreased. Thus, a heterogeneous subset should not receive so many labeled instances. In our case, a heterogeneous subset may be split and the risk improved, it is thus necessary to get more labeled instances in those subsets.

The value of the parameter and thus the true risk is unknown. As the selection strategy takes the form of a maximum of an unknown value, it falls in the scope of Optimism in the Face of Uncertainty. Under this paradigm, confidence intervals must be computed for each criterion, and the subset to sample is chosen according to the maximum upper bound of the intervals. In our problem, as we study classification, the family of the label distribution is known. We can therefore use Bayesian Credible Intervals for the true risk.

The family of the label distribution being Bernoulli, the family of conjugate prior is Beta distributions. By Bayesian inference:

∀ k ∈ [0, 1] P(µ k ≤ k ) = I k (T s,t μs,t + 1, T s,t (1 -μs,t ) + 1),
with T s,t the count of labeled instances in subset s and I . (α, β) the incomplete Beta function of parameter α and β.

Thus, the cumulative distribution function of the true risk is:

∀ k ∈ [0, 1], P(R s,t ≤ k ) =    I k W k (T s,t μs,t + 1, T s,t (1 -μs,t ) + 1) if [μ s,t ] = 0 1 -I 1-k W k (T s,t μs,t + 1, T s,t (1 -μs,t ) + 1) if [μ s,t ] = 1.
Suppose given δ 2 ∈ [0, 1], we compute the upper bound k (δ 2 ) of the Bayesian Credible Interval as the value of for which the probability of the true risk being higher is equal to δ 2 :

P(R s,t ≥ k (δ 2 )) = δ 2 ⇐⇒ k (δ 2 ) = W k × . . . (7)
I -1 (1 -δ 2 , T s,t μs,t + 1, T s,t (1 -μs,t ) + 1) if [μ s,t ] = 0 (1 -I -1 (δ 2 , T s,t μs,t + 1, T s,t (1 -μs,t ) + 1)) if [μ s,t ] = 1,
where I -1 (., α, β) is the inverse of the incomplete Beta function of parameter α and β.

The selection strategy is thus to sample the subset k t for which the upper bound of the Bayesian Credible Interval for the true risk is maximum:

k t = arg max k k (δ 2 ).
In this section, we introduced an Active Learning algorithm based on online classification trees. Considering an adaptive partitioning of the instance space, it achieves to allocate more effort to regions needing most. The online classification tree is well fitted to the Active Learning of an adaptive partition because it does not allow to remove splits. Indeed, it is of major importance in the Active Learning algorithm in order to get meaningful estimates of the parameter. Thus, at any time step, the distribution of labeled instances inside a subset must follow the distribution of unlabeled instance. Which is not the case if a split is removed and one of the children was favored over the other during the Active Learning process.

V. EXPERIMENTAL SETUP

In this section, we design experiments to evaluate the online classification tree algorithm introduced in the previous sections and its extension to the Active Learning problem. We design two kinds of experiments. The first on is a toy problem, with the true distributions parameters fixed by the user. The second is on real-world problems, where the number of available instances is limited. a) Toy problem: This provides a way for a better understanding of the behavior of the algorithm, but also comes with several advantages which we explain here. This problem consists in a predefined binary tree, to each leaf is given the parameter of a Bernoulli distribution and a size. The parameters and sizes are propagated back up in the tree, in order to know the parameter of all the none leaf nodes.

While learning the online classification tree on this problem, the only tests available are those of the predefined tree. Thus, only one test is considered for each node. The purpose of the algorithm is thus only to decide if the test must be installed or not, and not to chose which one to install. One limitation is that the most refined tree that can be outputted by the algorithm is the whole predefined tree, but this one is tried to be taken sufficiently refined.

At each time step, the next sample comes from a leaf of the predefined tree, randomly drawn from a set of possible leaves with a probability in proportion of its size. In the online learning scenario without Active Learning, the set of possible leaves always contain all the leaves of the predefined tree. In the Active Learning scenario, the set of possible leaves contains only the leaves descending from the subset selected by the Active Learning algorithm. The new label is then drawn from the corresponding Bernoulli distribution for this subset.

The use of known distribution parameters allows avoiding using a test set for evaluation. Indeed, the true risk is entirely known by the evaluation procedure (but not by the algorithm). Also, the number of samples that the algorithm can receive is not limited. Their is no need of instances, as the labels are given to subsets of the predefined tree, and the size of all subsets are given. Consequently, this problem is independent of the number of attributes. Indeed, it is only defined by the predefined tree, which can adapt to any instance space.

In our experiments, the toy problem has been fitted to an instance space which is [0, 1] 2 being split in two along alternating dimensions, until the number of leaves reaches 2 6 . The parameters of the leaves were set to

µ k = σ(5(x k,2 - √ x k,1 )),
where x k,1 and x k,2 represents the coordinates of the center of subset k, and ∀x ∈ R, σ(x) = 1 1+exp -x is the sigmoid function. The size of each leaf were equal. b) Real-world problems: In our experiments, we use several datasets from the UCI Machine Learning repository [START_REF] Lichman | UCI machine learning repository[END_REF]. Before each run of an experiment, the dataset is divided into two equal sets. The first is used as a pool from which the instances and their corresponding labels used for training are drawn. The unlabeled instances in the pool also serves to compute the size of each subset. The second is the test set, used for evaluation. At each time step, the algorithm makes predictions for instances in the test set, and the number of misclassified instances is recorded. The labels in the test set are totally hidden from the algorithm, which does not get any feedback from its prediction.

Experiments consisted into a certain number of runs of the algorithms. The performance in each run were averaged. The toy problem experiment consisted in 100 runs of the algorithms and a budget of 1,000 labeled instances. The realworld experiments consisted in 1,000 runs of the algorithms and was not stopped until all the instances in the training pool were labeled.

Our algorithms are compared to the following methods:

• c4.5 [START_REF] Quinlan | C4. 5: programs for machine learning[END_REF]: This is the most common batch learning classification tree algorithm. It uses information entropy as a measure of heterogeneity.

• Hoeffding: This is a modified version of VFDT [START_REF] Domingos | Mining high-speed data streams[END_REF] where the best split is not compared to the second best but only to no-split. Indeed, VFDT needs many instances before splitting in order to be sure that the installed test is the one that would be installed under perfect information. This is impractical for datasets with a limited number of instances. Hence, the only difference with BOCT is the use of Hoeffding bounds, which we want to compare.

• mbinomial [START_REF] Collet | Active learning for classification: An optimistic approach[END_REF]: This is an Active Learning algorithm which requires a fixed partition. The choice of the partition is non trivial for a real-world dataset. For this reason, we only test it on the toy problem. We instantiate this algorithm with two different partitions.

The first partition we use is whole predefined tree, we denominate the resulting algorithm by: mbinomial_64.

The second partition we use is a subtree of depth 3 of the predefined tree, we denominate the resulting algorithm by: mbinomial_8.

For all the algorithms used for comparison and for each datasets, the parameter δ (δ 1 and δ 2 in the case of BOCT-AL) has been tuned using a gridsearch.

VI. EVALUATION

In this section, we focus on comparing the algorithms introduced in the previous sections to state-of-the-art algorithms. For this, we use benchmarks described in section V. In order to evaluate the performance, we plot the true risk, averaged over the runs, function of the number of labeled examples. c) Toy problem: We start by showing off the results for the toy problem, results are displayed on Fig. 5. The goal of this experiment is to show the effect of using an adaptive partitioning of the instance space rather than a fixed partition. The mbinomial algorithm requires a fixed partition which cannot change during the process. The more refined the partition, the best optimal performance. We can see that with 64 subsets of equal size in the partition, mbinomial achieves a better performance at the end than with only 8 subsets. However, having a great number of subsets affects the learning rate. Indeed, in order to get an equivalent precision about the estimates of the parameters in each subsets, the number of points needed is in proportion of the number of subsets. This is confirmed by the fact that the performances with a few labeled instances is better for mbinomial with 8 subsets than with 64.

By using an adaptive partition, the algorithm is allowed to start with a small number of subsets, and thus get a good quality of prediction for each one of them with few samples. When more samples are received, the algorithm may use a more refined partition, while keeping the same quality of prediction. BOCT and BOCT-AL achieve to split at the right time, by splitting a subset as soon as the gain in true risk resulting from predicting different labels in each children subset counters the loss in quality of the predictions due to getting fewer labeled instances in each one of them. We see that, this way, BOCT-AL gets the same performance as the best Finally, we can observe that BOCT-AL shows better performance than BOCT. Thus, the Active Learning strategy used in BOCT-AL is effective. It manages to select subsets for which either the prediction can be improved or which may be split when acquiring more labeled instances.

To conclude on this experiment, we designed an adaptive partition algorithm which extends the mbinomial algorithm, by letting the partition to be automatically chosen by the algorithm. It achieves to use at any time the best partition, which leads to best performance. Now that the toy problem experiment demonstrated that our adaptive partition algorithms can effectively replace fixed partition ones, the goal is now to test it on real-world datasets. d) Real-world datasets: In this experiment, we use four different datasets coming from the UCI Machine Learning repository [START_REF] Lichman | UCI machine learning repository[END_REF]. The results for each of these datasets are displayed on Fig. 6[a-d]. The objective of this experiment is twofold. First, we want to test the utility of using Bayesian Credible Intervals rather than using much simpler and general bounds such as Hoeffding's bounds. Second, we want to measure the effect of actively selecting the labeled instances on the performances.

We can see that the algorithm using Hoeffding's bounds fail to achieve performances comparable to either BOCT or c4.5 on three out of four datasets. On the contrary, when comparing BOCT to c4.5, we observe that the final risks are equal on two datasets and higher for BOCT on two datasets. In addition, on Diabetes, where the Hoeffding's version is as good as c4.5, BOCT is better. Thus, BOCT is always better than the Hoeffding's version. Indeed, the theory on which the Hoeffding's version is built also applies to BOCT, but the Bayesian Credible Intervals are much tighter than Hoeffding's bounds. Hoeffding's bounds are known to be very tight when the number of samples is sufficiently large. However, for an online classification tree, we may want to split as soon as possible, which means when having received only few labeled instances. It is thus of major importance to used bounds that are tight with only few samples.

One other observation that we can have on these figures is that, even though BOCT attain final performances comparable to c4.5 and even better, the initial performances are always better for c4.5. The reason for this relies on the batch nature of c4.5. This way it does not have to plan ahead the splits of the instance space. At each time step, it considers the current best splits. For example, when having received only two labeled instances, one of each label, it may split the instance space to separate them, with a good chance that it increases the performance. On the contrary, in the online scenario, as no split nodes can be removed of the tree, an early bad split affects the quality of the prediction as it requires more labeled instances to get an equal accuracy. The algorithm may want to wait before splitting, as it may lead to a better split. During this waiting period, the performance is necessarily lower than in the batch scenario, but this enables a comparable long term performance.

With all the instances being labeled, we can see that the batch algorithm performs worse than an online algorithm. The reason for this is that in c4.5, the measure of heterogeneity is totally confident on the observations. Whereas in BOCT, an uncertainty is considered through the confidence intervals. Thus, c4.5 splits too much noisy subsets (parameter of Bernoulli close to 0.5) containing few labeled instances.

One may notice that the c4.5 curves are noisier than other algorithms' curves. In fact, this noise is also present on the other algorithms but only over the runs. In the batch scenario, at each time step, the classifier is recomputed, which makes the noise to appear on the curve. In the online learning scenario, the classifier is only updated, thus, the noise is the same all along one specific run.

We see that BOCT-AL shows slightly better performance than BOCT. We can notice that BOCT-AL has exactly the same performance as BOCT at the beginning of the process. This is because before the first test has been installed, there is only one subset, and the Active Learning is useless since it serves to select one subset. We also observe that at the end, the performances of BOCT-AL and BOCT are the same. Indeed, even though the priority of one cluster among another is different, the decisions in a particular subset are not affected by the state of other subsets. This means that Active Learning only changes the order in which subsets are studied. When all the instances in the pool have been labeled, the Active Learning does not have an effect.

To conclude on this experiment, we have seen that using Bayesian Confidence Intervals instead of Hoeffding's bounds greatly improves the performance of the algorithm. In the case of a large number of instances to label, when anticipating on the future is necessary, an online classification tree is preferred to a batch one. And the BOCT performances can be further improved by using the Active Learning criterion designed in this paper.

VII. CONCLUSION

In this paper, we show that Bayesian Credible Intervals are well suited for learning an online classification tree, as the family of the output distribution is known. This enables using confidence intervals which are perfectly tight, both with few and many samples. We introduce BOCT and BOCT-AL, two novel algorithms for Online and Active Learning of a classification tree. The first one uses Bayesian Credible Intervals on the heterogeneity of labels to decide when to split. The second one adds the use of Bayesian Credible Intervals on the risk in order to select the instances to be labeled. Both the algorithms are evaluate on a toy problem, as well as on several real-world datasets, and compares positively to stateof-the-art. Following work may focus on a theoretical analysis of these algorithms, this could lead to an automatic choice of the parameter. This works focus on binary classification, this could be easily adapted to the multi-class problem. Finally, such trees could be used in a forest of trees [START_REF] Breiman | Random forests[END_REF].

ACKNOWLEDGMENT

Figure 1 :

 1 Figure 1: BCI on weighted sum of true heterogeneities after several tests compared to one of the parent node, here, test 4 would be installed.

Figure 2 :

 2 Figure 2: Core algorithm

Figure 3 :

 3 Figure 3: UpdateTree

Figure 4 :

 4 Figure 4: Active Learning algorithm

Figure 5 :

 5 Figure 5: Toy problem

  The InterCellMISN program of the French State to Lorraine region 2007-2013 aims at setting up a cluster (256 PCs) for interactive fine grain computation. It is granted by INRIA and the Region Lorraine (CPER 2007). The experiments reported in this paper where lead on this hardware.

  Figure 6


F Orieux
email: orieux@pasteur.fr.

J.-C Olivo-Marin

E Sepulveda

V Loriette
email: vincent.loriette@espci.fr.

B Dubertret
email: benoit.dubertret@espci.fr.

Bayesian Estimation for Optimized Structured Illumination Microscopy

Franc ¸ois Orieux * , Eduardo Sepulveda, Vincent Loriette, Benoît Dubertret, Jean-Christophe Olivo-Marin * Abstract-Structured Illumination Microscopy is a recent imaging technique that aims at going beyond the classical optical resolution by reconstructing high resolution images from low resolution images acquired through modulation of the transfer function of the microscope. The classical implementation [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF] has a number of drawbacks, like requiring a large number of images to be acquired and parameters to be set manually in an ad-hoc manner that have, until now, hampered its wide dissemination. Here, we present a new framework based on a Bayesian inverse problem formulation approach that enables the computation of one high-resolution image from a reduced number of low-resolution images and has no specific constraints on the modulation. Moreover, it permits to automatically estimate the optimal reconstruction hyper-parameters and to compute an uncertainty bound on the estimated values. We demonstrate through numerical evaluations on simulated data and examples on real microscopy data that our approach represents a decisive advance for a wider use of high-resolution microscopy through structured illumination.

Index Terms-Structured illumination microscopy, image reconstruction, microscopy, deconvolution, super-resolution, Bayesian inversion, MCMC.

I. INTRODUCTION

A. Superresolution Microscopy

O PTICAL microscopy is one of the major research tools of modern biology as it opens the way to study the links between cellular phenotypes and the underlying molecular mechanisms through the quantitative imaging of cell compartment dynamics [START_REF] Galy | Nuclear pore complexes in the organization of silent telomeric chromatin[END_REF], [START_REF] Arhel | Quantitative four-dimensional tracking of cytoplasmic and nuclear hiv-1 complexes[END_REF]. Classical optical microscopes are, however, limited to a lateral resolution of 200 nanometers and a 600 nanometer axial resolution imposed, by the phenomenon of diffraction of light by the microscope objective pupil, which does not allow for resolving a large class of sub-resolution biological objects.

Two active fields of research aim at advancing the capabilities of optical microscopes either by increasing the resolution or by improving the localization accuracy of individual objects [START_REF] Vonesch | The colored revolution of bioimaging[END_REF], [START_REF] Hell | Microscopy and its focal switch[END_REF]. Examples of the latter are recent techniques like PALM [START_REF] Betzig | Near-field Scanning Optical Microscopy[END_REF] or STORM [START_REF] Huang | Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy[END_REF]: individual objects are excited with a high probability of being physically separated from their neighbors and therefore, a highly accurate estimation of their position can be computed. An image is constructed by accumulating individual objects' positions. It should be stressed that strictly speaking, no fundamental resolution enhancement is achieved in this case as it is not possible to distinguish two objects closer than the classical diffraction limit had they been excited simultaneously.

In the former approaches, the idea is to modify the sample illumination process in such a way that the classical diffraction limit does not hold anymore and higher spatial frequencies can be recovered, which means a true gain in resolution. For example, confocal microscopy [START_REF] Brakenhoff | Visualization and analysis techniques for three dimensional information acquired by confocal microscopy[END_REF], besides enabling 3D imaging of biological samples, has improved the lateral resolution by a factor of 1.4 by using conjugated pinholes to reject out-offocus light. Similarly, two-photon microscopy [START_REF] Denk | Two-photon laser scanning fluorescence microscopy[END_REF] and more recently STED [START_REF] Hell | Three-photon excitation in fluorescence microscopy[END_REF] have also increased the resolution by restricting the excitation volume. As powerful as they are, these techniques present the drawback of making inefficient use of photons as the ratio between the number of detected fluorescence photons to the number of excitation photons is very low. As an alternative to these scanning approaches, full field techniques like 4-pi microscopy [START_REF] Hell | Three-photon excitation in fluorescence microscopy[END_REF] and Structured Illumination Microscopy (SIM) [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF] have better photon efficiency because most of the photons are collected and used to form the image. Even if the resolution power of SIM approach is lower compared to STED, the good efficiency and the wide field illumination constitute an important counterpart.

B. A signal processing interpretation of SIM

Diffraction theory states that incoherent1 optical systems can be described by the impulse response (or point spread function -PSF) obtained from the square of the Fourier transform of the aperture [START_REF] Goodman | Introduction to Fourier Optics[END_REF]. The corresponding transfer function (OTF) is the auto-correlation of the aperture and is equal to zero for all frequencies beyond a cutoff frequency f c . In this case, all the information outside this bound is completely lost. This phenomenon is illustrated in Fig. 1 where all the frequencies outside the support (in dotted black) of the transfer function are not observed.

The idea of structured illumination microscopy (SIM) [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF], [START_REF] Gustafsson | Threedimensional resolution doubling in wide-field fluorescence microscopy by structured illumination[END_REF], [START_REF] Mandula | Patterned excitation microsopy[END_REF] is to inject high-frequencies inside the support of the transfer function with the help of amplitude modulation before convolution. In other words, the objective of the illumination is to introduce aliasing through modulation, as illustrated in Fig. 1, and to gain information in the frequency range beyond f c in order to reconstruct a high resolution (HR) image.

Let us denote the original image f (x, y) ∈ L 2 where (x, y) ∈ R 2 and its continuous Fourier transform F (f x , f y ) (also called spectrum). In the standard SIM setup, the illumination modulation pattern is

m(x, y) = I 0 1 + α cos 2π(k x x + k y y + φ) (1)
and takes its values in [I 0 (1α), I 0 (1 + α)] (the illumination, being an intensity, is necessarily positive). The parameter α is the modulation depth, whose value depends on the coherent or incoherent nature of the emission light, on the aberrations induced by the observed sample and the overall quality of the optical setup. Its Fourier transform is

M (f x , f y ) = I 0 (δ(f x, f y )+ α 2 δ(f x -k x , f y -k y )e -2iπφ + α 2 δ(f x + k x , f y + k y )e -2iπφ . (2)
With this modulation, the main spectrum F (f x , f y ) is replicated three times and centered in (0, 0), (k x , k y ) and (-k x , -k y ). Let us now introduce the convolution with the PSF h(x, y), corresponding to a window in the Fourier domain with the transfer function H(f x , f y ) as

G(f x , f y ) = I 0 H(f x , f y ) (F (f x , f y )+ α 2 F (f x -k x , f y -k y )e -2iπφ + α 2 F (f x + k x , f y + k y )e -2iπφ (3)
where G(f x , f y ) represents the data in the Fourier domain.

Since the transfer function H is not replicated, the weighting coefficients are different for each term in [START_REF] Arhel | Quantitative four-dimensional tracking of cytoplasmic and nuclear hiv-1 complexes[END_REF]. This concept, illustrated for one specific orientation in Fig. 2a, allows to measure frequencies present in an extended support. By repeating the same modulation, but with different orientations, a tiling can be constructed giving raise to measurements in a full ring outside the classical support of the PSF as illustrated in Fig. 2b.

f f m -f c D - D + D 0

C. Existing approaches

The mere fact of introducing aliasing as done in Eqn. [START_REF] Arhel | Quantitative four-dimensional tracking of cytoplasmic and nuclear hiv-1 complexes[END_REF] is not sufficient to obtain a HR (High Resolution) image. Indeed, with just one LR image (Low Resolution), aliased frequencies are summed up and there is no way to solve the

f x f y (a) f x f y (b)
Fig. 2. Measure of high-frequency by modulation. Fig. 2a illustrates the modulation with the reproduction, in gray, of the central pattern around the two Dirac (black dot). Fig. 2b illustrates the tiling for the measure of highfrequencies in a second ring. ambiguity problem unless making use of supplementary information. This can be done either by acquiring additional LR images with different illumination patterns, or by a different modulation.

In the original method [START_REF] Gustafsson | Threedimensional resolution doubling in wide-field fluorescence microscopy by structured illumination[END_REF]- [START_REF] Conchello | Image estimation for structured-illumination microscopy[END_REF], each pattern centered on (0, 0), (k x , k y ) and (-k x , -k y ) is considered to be a different unknown. From Eqn. (3) we can write

G(f x , f y ) = I 0 (H(f x , f y )F (f x , f y )+ α 2 H(f x , f y )F (f x -k x , f y -k y )e -2iπφ + α 2 H(f x , f y )F (f x + k x , f y + k y )e -2iπφ (4)
= H(f x , f y ) (F (f x , f y ) + F -(f x , f y ) + F + (f x , f y )) . (5) = D 0 (f x , f y ) + D -(f x , f y ) + D + (f x , f y ). (6)
Aliasing appears through the sum of the three components D 0 , D -and D + , and in this case, the resolution of the system consists in separating the three components for each orientation. By observing that the angular term -2iπφ appears only in the replication D -and D + , the classical solution relies on using three different patterns for each orientation, which are shifted in phases relative to one another by 2π/3. We can, therefore, rewrite (6) as:

G k (f x , f y ) = D 0 (f x , f y ) + D -(f x , f y )e -2kiπ/3 + D + (f x , f y )e 2kiπ/3 . ( 7
)
where k = 0, 1 or 2. Thereafter, thanks to the choice of the 2π/3 phase shift, the system can be solved by linear combinations between the LR images as:

D 0 = 1 3 (G 0 + G 1 + G 2 ) (8)
D -= 1 3 G 0 + e 4iπ/3 G 1 + e -4iπ/3 G 2 (9)
D + = 1 3 G 0 + e -4iπ/3 G 1 + e 4iπ/3 G 2 . ( 10
)
Following this component separation for one orientation, the classical approaches by Gustafsson [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF] or Mandula [START_REF] Mandula | Patterned excitation microsopy[END_REF] go through a number of different processing steps:

• registration of sub-patterns,

• correction of the weighting by transfer function,

• new combination of the sub-patterns with the new weighting factors, • final reconstruction through Wiener filtering.

Overall, these methods suffers from several drawbacks. First, the use of the Wiener filter is problematic: to regularize the problem, it does require prior information, which usually involves at least an SNR (Signal to Noise Ration) parameter unknown in real-life problems. Effectively, the SNR parameter depends on the noise power, that can be estimated, but also on the unknown signal. Consequently, this regularization parameter is generally tuned in an ad-hoc way. Also, this approach requires at least three images per modulation orientation, with a total of nine LR images to reconstruct a single HR image. Finally, the existing approaches do not provide a rigorous formalism with which to derive an uncertainty bound on the estimation.

D. Development of this paper

In this paper, we propose a new global and coherent Bayesian framework, to compute the high-resolution (HR) image 2 , which is able to address all the previous drawbacks. It relies on an inverse problem formulation approach that enables the computation of one HR image from just four low-resolution images and without any constraint on the modulation pattern (like for example, the 2π/3 phase shifting factor). Moreover, within this Bayesian framework [START_REF] Demoment | Image reconstruction and restoration: Overview of common estimation structure and problems[END_REF], it is possible to rigorously formulate prior information about the system and to estimate jointly all the unknown regularization parameters. Finally, the Bayesian framework also permits to estimate an uncertainty bound about the reconstructed image and the associated parameters. The proposed method is more objective and requires less or no tuning by the user.

The paper is organized as follows. The proposed models and algorithms are presented in Section II. Section III discusses aspects related to the number of LR images. Section IV is devoted to experiments and results.

II. BAYESIAN INVERSION FOR SIM

The reconstruction of an HR image can be seen as an inverse problem, i.e., the estimation of the original image that is at the origin of the data. There is abundant literature about this subject [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF]- [START_REF] Demoment | Inverse problems: theory and applications[END_REF].

A powerful framework to address the inversion issue is the Bayesian approach. In this framework, all the available information about the unknowns [f , γ] in an uncertain environment is contained in the posterior law p(f , γ|g), where f is the variable of interest, γ the hyper-parameters and g the data. By applying the Bayes rules, we get p(f , γ|g) = p(g|f , γ)p(f , γ)/p(g). [START_REF] Goodman | Introduction to Fourier Optics[END_REF] This posterior law can be obtained thanks to the product of the likelihood p(g|f , γ), that models all the information coming from the data and their uncertainty, and the prior law of the unknowns p(f , γ) that models the information coming from other sources (like the classical smoothness). This is an inverse formula that allows determining the unknown causes (the image f ) from the known effects (the data g).

A. Forward model

The forward model aims to simulate or reproduce the L LR data images acquisition g l ∈ R N , indexed by l, from the original HR image f ∈ R N . The original image f is written as a vector where the lines of the image are column stacked.

As described in Sec. I, there are two main elements to be modeled: the illumination and the convolution by the optical setup.

The illumination is an amplitude modulation by an excitation light. Each pixel of the image or each element of f is weighted by a cosine pattern. This operation is described by the multiplication of f with a diagonal matrix M l where the diagonal elements correspond to the modulation pattern of Eq. ( 1). This matrix depends on the modulation parameters, different for each l LR image. Since we consider the case of an ideal illumination model Eq. ( 1), the influence of the optical setup reduces to a change in the contrast of the grid, easily modeled by the α parameter described previously and whose value is fixed by the optical setup and the sample. For coherent illumination, the OTF is flat and has therefore no influence on α. For incoherent illumination, the OTF decreases with frequency and the value of α is attenuated, whereby fringes become less visible. In any case, α must be either fixed in an ad-hoc manner or estimated from the data.

After modulation, the optical system collects the light and focuses it on a detector. The lens and the detector global response can be modeled, in the incoherent light regime, as a convolution that combines both responses. The convolution operation is written as a product between the modulated image and the matrix H. This is a Toeplitz matrix and its first line corresponds to the discrete 2D PSF H.

These two previous operations are combined to model the full acquisition. For the acquisition of one data image we have

g l = HM l f + n l ( 12
)
where n l is an unknown term that models the model error and the noise.

If the matrix H is block-circulant circulant-block (BCCB), it can be diagonalized in Fourier domain F † Λ H F = H, where F is the unitary Fourier transform and Λ H a diagonal matrix. This is the discrete formulation of the duality of the convolution in direct and Fourier domain. In this case, the model becomes

g l = F † Λ H F M l f + n l . ( 13
)
By collecting all the L LR images in one vector g = [g 1 , g 2 , . . . g L ] (same thing for n), the forward model can be written as

g = HM Rf + n ( 14
)
where R = [I t I t . . . ] t is a replication matrix, M a block diagonal matrix with M l in each block and H a block diagonal matrix with H in each block.

An important issue about this forward model is that the matrix M can model any kind of amplitude modulation, allowing a very flexible framework to design efficient illumination. In particular, this modeling allows us to reduce the number of necessary images as described in the next section and illustrated in Sec. IV.

B. Likelihood

The noise term n is known up to statistical information. The choice of the noise law depends mainly on the application and most cases in image processing use either Gaussian or Poisson distributions. In the case of position dependent noise, Poisson or non-homogeneous Gaussian distributions can be used with the difficulty that in the former case the possible estimators are no longer linear wrt. data, while in the latter the covariance matrix needs to be estimated through quite cumbersome procedures.

In applications where a high number of photons are collected, some simplifications occur. First, the Poisson law tends towards a Gaussian. Then, the noise introduced by readout electronics can be described accurately by a Gaussian statistics. Finally, even if the standard deviation depends on the position, a Variance Stabilizing Transform like in [START_REF] Zhang | Wavelets, ridgelets and curvelets for poisson noise removal[END_REF] can be used to stabilize it. For these reasons, our proposed method uses an iid Gaussian law

n l ∼ N 0, γ n -1 l I ( 15
)
where γ nl is the unknown inverse variance for this LR image. Consequently, the data law when the HR image is known, or the likelihood, is also Gaussian:

p(g l |f , γ nl ) = (2π) -N/2 γ n N/2 l exp -γ nl g l -HM l f 2 /2 (16) = (2π) -N/2 γ n N/2 l exp -γ nl • g l -Λ H F M l f 2 /2 (17)
thanks to the Parseval relationship, with

• g = F g. The law for all the images is p(g|f , γ n ) = (2π) -N L/2 |Σ| -1/2 exp - 1 2 g -HM Rf t Σ -1 n g -HM Rf (18)
where Σ -1 n is block diagonal and each block is γ nk I. If the noise power is considered to have the same value for each image, then Σ -1 n = γ n I and the data law is

p(g|f , γ n ) = (2π) -N L/2 γ n N L/2 exp -γ n g -HM Rf 2 /2 . ( 19
)

C. Prior law of the HR image

The definition of the prior law for the HR image gives an opportunity to introduce additional information, e.g. a positivity constraint (not used in this paper). Since our forward model includes a convolution operator, some information on the smoothness features of the image must be introduced to cancel the amplification of noise, and to that goal, we use Gaussian fields.

The probability law for the HR image with a given covariance matrix Σ γ parametrized by a vector γ reads

p(f |γ) = (2π) -N/2 |Σ γ | -1/2 exp - 1 2 f t Σ γ f . ( 20
)
For computational efficiency, the precision matrix is designed (or approximated) in a toroidal manner, and is diagonal in the Fourier domain Λ f = F Σ γ F † . Thus, the law for f writes

p(f |γ f ) = (2π) -N/2 det[F ] det[Λ f ] 1/2 det[F † ] exp - 1 2 f t F † Λ f F f (21) = (2π) -N/2 det[Λ f ] 1/2 exp - 1 2 • f † Λ f • f . ( 22
)
This approximation is sometimes referred to as a Whittle approximation [START_REF] Calder | Introduction to Whittle (1953) 'The analysis of multiple stationary time series[END_REF] (see also [START_REF] Orieux | Bayesian estimation of regularization and PSF parameters for Wiener-Hunt deconvolution[END_REF] or [START_REF] Molina | Blind deconvolution using a variational approach to parameter, image, and blur estimation[END_REF]) for the Gaussian law.

In this paper, we focus mostly on smooth images, that can be modeled by positive correlations between pixels. This can be introduced by high-frequency penalties using any circulant differential operator (like Laplacian). The differential operator is denoted by D and its diagonalized form by Λ D = F DF † . Then, the inverse variance matrix writes Σ -1 γ f = γ f D t D and its Fourier counterpart writes Λ f = γ f Λ D † Λ D where γ f is a positive scale factor tuning the degree of smoothness. The Wiener filter uses the identity matrix I instead of D.

In addition to the efficiency of the numerical computation due to Gaussian hypothesis, this model also offers as an advantage the knowledge of the normalization coefficient and its dependency on the hyper-parameter γ f . With this form of correlation, the law can be written as:

p(f |γ f ) ∝ γ f (N -1)/2 exp - γ f 2 Λ D • f 2 . ( 23
)
The parameters γ n and γ f are usually tuned empirically. The next section presents a way to determine them automatically.

D. Prior law of hyper parameters

The parameters γ n and γ f are considered unknown. With respect to the Bayesian framework

p(f , γ n , γ f |g) ∝ p(f |g, γ n )p(f |γ f )p(γ n )p(γ f ) (24)
prior laws p(γ n ) and p(γ f ) for the parameters must be defined. A classical choice for hyper parameter laws relies on conjugate prior [START_REF] Robert | The Bayesian Choice[END_REF]: the conditional posterior for the hyperparameters is in the same family as its prior. All the parameters are precision parameters of a Gaussian law and a conjugate law for these parameters is the Gamma law. Given parameters

(α i , β i ) the pdf reads p(γ i ) = G(α i , β i ) (25) = 1 β αi i Γ(α i ) γ αi-1 i exp (-γ i /β i ) , ∀γ i ∈ [0, +∞[. (26)
Very little prior information is available for these parameters.

With parameters set to (0, +∞) values, one obtains the Jeffreys' non-informative prior. Jeffreys' law is a classical law for the precisions and is considered as non-informative [START_REF] Kass | The selection of prior distributions by formal rules[END_REF].

E. Posterior mean estimator

The joint posterior law is defined as

p(f , γ n , γ f |g) ∝ γ f (N -1)/2-1 L-1 l=0 γ n N/2-1 l exp - L-1 l=0 γ nl 2 g l -HM l f 2 exp - γ f 2 Df 2 . ( 27
)
This law is multidimensional (the number of dimensions is the number of pixels in the HR image along with the number of unknown hyper-parameters) and quite complex. Clearly, it is not possible to manipulate or keep all the information embedded in the posterior law. The classical solution relies on the choice of a particular point like the point of maximum probability called maximum a posteriori (MAP)

f , γ n , γ f = arg max f ,γ n ,γ f p(f , γ n , γ f |g) (28)
that leads to an optimization problem, or the mean of the law

f , γ n , γ f = [f , γ n , γ f ] p(f , γ n , γ f |g) df dγ n dγ f (29
) that leads to an integration problem. Both estimators are possible, but the posterior mean has some advantages. First, the mean is known to minimize the mean square error (MSE) [START_REF] Robert | The Bayesian Choice[END_REF]. Second, use of classical algorithms to compute the mean allows access to all the other moments of the law, including the uncertainty of the estimate. For these reasons, the mean is computed with the help of a Gibbs algorithm.

F. Gibbs sampler

To compute the posterior mean, a Monte Carlo Markov chain is used to provide stochastic samples. The provided samples allow computing moments of the law, like the mean and the standard deviation, assimilated to uncertainty on the mean. Afterwards, samples are also used to approximate marginal laws as histograms.

The Gibbs algorithm [START_REF] Robert | The Bayesian Choice[END_REF], [START_REF] Winkler | Image Analysis, Random Fields and Markov Chain Monte Carlo Methods[END_REF], [START_REF] Brémaud | Markov Chains. Gibbs fields, Monte Carlo Simulation, and Queues[END_REF], first introduced in image processing by Geman & Geman [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF], is an MCMC algorithm to sample any law. It relies on iterative sampling of the conditional posterior law for a set of parameters given the other obtained through previous iterations. The structure of the algorithm is a repetition of following steps:

1) Draw f ∼ p f |γ n (k) , γ f (k)
, g , which is a multidimensional Gaussian law. The size and complexity of the covariance law prevents direct sampling. Classical algorithms like per-pixel sampling [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF] or Cholesky factorisation [START_REF] Rue | Fast sampling of Gaussian Markov random fields[END_REF], [START_REF] Lalanne | Stochastic artificial retinas: algorithm, optoelectronic circuits, and implementation[END_REF] are impossible in practice because of either the required time or of memory. The appendix A proposes an adapted algorithm of [START_REF] Rue | Fast sampling of Gaussian Markov random fields[END_REF], [START_REF] Lalanne | Stochastic artificial retinas: algorithm, optoelectronic circuits, and implementation[END_REF] approaches based on the resolution of an optimization problem.

2) Draw γ nl ∼ p γ nl |f (k+1) , γ f (k) , g , which is a Gamma law k+1) , g , also a Gamma law

p γ nl |f (k+1) , g = G N L/2, 2/ g l -HM l f (k+1) 2 , 3) Draw γ f ∼ p γ f |f (k+1) , γ n (
p γ f |f (k+1) = G (N -1)/2, 2/ Df (k+1) 2
.

After a burn-in or convergence time, the chain becomes stationary and the complete set of samples {f (k) , γ f (k) , γ n (k) } k∈K follows the joint posterior law, Eqn. [START_REF] Brémaud | Markov Chains. Gibbs fields, Monte Carlo Simulation, and Queues[END_REF]. Then, the posterior mean is approximated by an average of the samples

f = E[f ] ≈ 1 K K-1 k=0 f (k) (30)
using the large numbers law.

In addition to the mean, the samples can be used to compute other moments of the posterior law. For example, the posterior variance:

V[f ] = E (f -E[f ]) (f -E[f ]) t (31)
is used to provide uncertainty about the estimated value. In the case of the image, the covariance matrix is very large, but the diagonal value (variance of each pixel) is obtained with

V[f ] ≈ 1 K K-1 k=0 f (k) f (k) -f 2 (32)
where is the term-wise multiplication.

III. NUMBER OF LR IMAGES

In the standard SIM approach [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF], [START_REF] Mandula | Patterned excitation microsopy[END_REF], each sub-pattern in Fig. 1 is considered to be a different unknown. The particular illumination structure with 2π/3 phase shifting allows the resolution of the system simply by linear combination between nine images. The idea to build an equation system that takes into account the redundancy in the data and achieves a reduction of the number of required images is presented in [START_REF] Heintzmann | Saturated patterned excitation microscopy with twodimensional excitation patterns[END_REF]. In that paper, the possibility to reduce the number of images down to four images in linear SIM is already stated and thought of, however without specific implementation details and examples.

Here, thanks to an inverse problem approach, we prove an effective realisation of the use of redundancy and show that in some conditions and depending on the optical setup, the number of images can indeed be decreased to as low as four.

The rationale is based on several points. First, in the standard approach, the D 0 sub-pattern is reconstructed three times, once for each orientation, although it is the same image in all three cases 3 In our proposed approach,it is not necessary to acquire additional images, because the presence of D 0 is naturally taken into account by the joint inversion of all the data. This effectively allows to reduce the number of image to 7.

Second, the original image is real, and consequently the pattern F + and F -, in Eq. ( 5), are identical up to a conjugacy. Moreover, the satellite and central patterns are equal on some interval F + (f mf ) = F * (f ). Consequently, as illustrated in Fig. 3, as long as f m ≤ f c , there are only two unknowns inside the support of the OTF :

F (f ), ∀f ∈ [0, f c ] and F (f ), ∀f ∈ [f c , f c + f m ],
if and only if f m ≤ f c . In Fig. 3 this illustrated by U 0 and U 1 being identical. The extension to 2D is straightforward. Eqn. ( 6) can be rewritten as

f f m = f c U 0 U 2 U 1

OTF support

G(f x , f y ) = D 0 (f x , f y )+D + (f x , f y )+D * + (f x -2k x , f y -2k y ) (33)
and for each orientation, there are only two orders, D 0 and D + , instead of three. For one orientation we have exactly two unknowns if k 2

x + k 2 y ≤ f c , as illustrated in Fig. 4a, except for the surfaces in gray where three unknowns are present. When tiling the plane with images acquired with different orientations, some frequencies can be resolved, as illustrated in Fig. 4b. Frequencies where further ambiguities remain can be estimated thereafter with the information coming both from the likelihood (which is not sufficient by itself ) and the prior law.

To solve the image reconstruction problem with a lower data redundancy, the approach is to have an equation system that considers all the data jointly to estimate all the unknowns, without ambiguity, as discussed in [START_REF] Heintzmann | Saturated patterned excitation microscopy with twodimensional excitation patterns[END_REF]. Our proposed inversion approach, which uses a direct model and all the data jointly to reconstruct the image in one global step, is an implementation of the latter. If we consider, for the sake of simplicity, the same noise level for all image, then the step 1) 3 This is true when a sufficiently high density of photons with random polarisation are acquired, leading to an indistinguishable polarisation. This was the case for our data: we had a minimum of 450 photons per pixel, and the conditions were met. It should be noted that birefringent samples would however be a limitation of the setup. in the algorithm (Sec. II-F) corresponds to the simulation of the conditional posterior law:

p(f |g, γ n , γ f ) ∝ exp - γ n 2 g -HM f 2 - γ f 2 Df 2
The Gaussian law can be handled through a quadratic criterion (see Appendix. A)

J(f ) ∝ γ n g -HM f 2 + γ f Df 2 .
where the optimisation leads to the resolution of the equation system

M HH t M t + D t D f = M t H t g (34)
effectively solved with a conjugated gradient descent algorithm (see [START_REF] Nocedal | Numerical Optimization, ser. Series in Operations Research[END_REF] and Appendix. A). In the equation system Eq. [START_REF] Rodríguez | Axial coding in full-field microscopy using three-dimensional structured illumination implemented with no moving parts[END_REF] all the pixels of the image f are considered jointly with all the data collected in g.

IV. EXPERIMENTS

In this section, the proposed method for super-resolution imaging with structured illumination is studied. The study focuses mainly on image quality and on the contributions of the proposed approach.

A. Simulations

The experiments are conducted on two simulated cases, "Barbara" and the "test pattern", well suited to analyze the improvement in resolution and on real data obtained with "Structured Illumination Microscopy" [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF], [START_REF] Gustafsson | Threedimensional resolution doubling in wide-field fluorescence microscopy by structured illumination[END_REF], [START_REF] Mandula | Patterned excitation microsopy[END_REF]. The results are also compared with a classical deconvolution approach [START_REF] Mandula | Patterned excitation microsopy[END_REF], [START_REF] Orieux | Bayesian estimation of regularization and PSF parameters for Wiener-Hunt deconvolution[END_REF].

The data has been simulated by our forward model described in Sec. II-A. Without loss of generality, the PSF is the classical Airy disc often encountered in optical problems [START_REF] Goodman | Introduction to Fourier Optics[END_REF]. Its expression in Fourier space is

H(f ) = 1 π 2 cos -1 f f c -sin 2 cos -1 f f c , (35)
∀|f | ∈ [0, f c ], 0 otherwise. With this PSF, the cutoff frequency is at f c , beyond which no information is available.

The modulation patterns are described in Sec. I-B Eq. ( 1). The number of LR images is four with three different orientations to construct a "second ring" of frequency. The orientations are -30, 30 and 90 degrees. In all cases, the reconstruction of the HR image is done with only four LR images: one without structured illumination and three with illumination in the three different orientations.

The regularity of the reconstructed image is measured through the Laplacian D = [0, -1, 0; -1, 4, -1; 0, -1, 0] and Λ D is obtained with an FFT-2D. Regarding the parameters of the Gamma laws, we set α and β to (0, +∞) in order to choose the non-informative Jeffreys' law [START_REF] Kass | The selection of prior distributions by formal rules[END_REF].

The algorithm used is described in Sec. II-F. A burn-in period of 20 samples is set before the computation of the empirical mean by Eqn. [START_REF] Lalanne | Stochastic artificial retinas: algorithm, optoelectronic circuits, and implementation[END_REF]. The algorithm 1 is stopped following the criterion

f (K+1) -f (K) 2 / f (K) 2 ≤ , (36)
when the difference between two successive means is less than = 10 -3 , where K is total number of samples.

The results are compared to the classical approach [13] that requires nine LR images. To observe the improvement in image quality provided by structured illumination, the results are also compared to a multichannel unsupervised deconvolution approach like [START_REF] Orieux | Bayesian estimation of regularization and PSF parameters for Wiener-Hunt deconvolution[END_REF] with four channels (4 LR images). Since the noise model is white Gaussian, this corresponds to a deconvolution with 1 LR image and 4γ n higher precision level.

The objective metrics used to quantify the quality of the results are (i) the improved signal-to-noise ratio (ISNR) of estimate f 1 over f 2 defined as

ISNR f 1 , f 2 = 10 log 10 f -f 2 2 / f -f 1 2 (37)
and (ii) the SSIM (Structural SIMilarity) [START_REF] Wang | Mean square error : Love it or leave it ? A new look at signal fidelity measures[END_REF] designed to be adapted to human visual perception

SSIM f 1 , f 2 = (2µ 1 µ 2 + c 1 )(2cov 12 + c 2 ) (µ 2 1 + µ 2 2 + c 1 )(σ 2 1 + σ 2 y + c 2 ) (38)
where µ and σ are the mean and standard deviation of the vectors, cov 12 the covariance between f 1 and f 2 and c 1 and c 2 are stabilisation terms. The SSIM index is a measure between an image and a reference (ground truth here) and takes its value between -1 and 1 (1 meaning full identity with the reference).

The algorithm has been applied on two sets of simulated data:

• The "Barbara" image, Fig. 5a, which is a 519 × 519 pixel image. The cutoff frequency is set to 0.3 in reduced frequency, and therefore almost all the frequency plane of the original sampled image can be observed. In this simulated image, frequencies above 0.5 are not present and therefore no expansion of the frequencies plane is required. • The "test pattern" image, Fig. 5b, which is a 263 × 263 pixel image. The cutoff frequency is set to 0.12 in reduced frequency. The tests have been conducted with the following conditions:

• SNR ≈ 80 dB. This is a favorable case, with standard deviation σ ≈ 1 for an intensity range in [0, 255], and approximately 6500 photons per pixel (for real data acquired with an ANDOR iXon 885 camera and a shotnoise model), and • SNR ≈ 20 dB. This is an unfavorable case, with σ ≈ 6 and a number of photons of approximately 450 per pixel, and where the same noise level is applied to all LR images. The real data present an SNR in the range 20 to 40 dB.

1) Estimation results: Results on "Barbara" are illustrated in Fig. 6. For comparison, the Fig. 6a is the data obtained with simulated brightfield illumination (without structured illumination) and 6e is a slice. They are the results of a simulated convolution with an optical system. Compared to the true image 5a, small details are no longer visible, especially, the stripes on clothes. Results from deconvolution are illustrated in 6b and 6f. Some details are restored around pixel 150, but the stripes are still not visible. Results from the structured illumination with the proposed approach are illustrated in Figs. 6c and6g. With this approach, many more small details are visible. The stripes on clothes are visible and high frequencies around pixel 25 or 225 are clearly restored in comparison to the deconvolution. Our results are virtually identical to the classical approach [START_REF] Mandula | Patterned excitation microsopy[END_REF], illustrated Figs. 6d and6h, where 9 images are necessary.

This gain in high frequencies is clearly visible in the "spectrum" illustrated in Fig. 7. Since the optical transfer function is zero outside f c , the deconvolution process is unable to restore frequencies outside this limit, as shown by Fig. 7c. The structured illumination, through modulation, is able to recover high frequencies above f c . This is well illustrated by Fig. 7d. The tiling pattern in the three orientations is visible. A close resemblance with the true spectrum in Fig. 7a can be noticed, as well as the presence of the high-frequencies necessary for reconstructing the stripes on clothes.

The fact that high frequencies are measured beyond the cutoff frequency f c of the optical transfer function is illustrated in the empirical estimation of the circular mean of the spectrum in Fig. 8a. The gain in low frequency with deconvolution is visible only up to f c , while It can be seen that structured illumination allows restoring frequencies well beyond that limit. In this experiment, the level of noise is low enough to allow a high quality restoration up to frequencies close to the highest one. Results with the "test pattern" in Fig. 9 are similar, notwithstanding lower modulation and cutoff frequencies, meaning that less high frequencies are measured. In comparison to the deconvolution result in Fig. 9b, the structured illumination in Fig. 9c increases the spatial resolution and visibility of small details, which is also obvious from the spectrum (Fig. 10) with the tiling. On the circular mean (Fig. 8b), it is clear that frequencies above f c ≈ 0.11 and up to f c + f m /2 ≈ 0.18 are restored. Numerical performances are reported in Tab. I. For a noise level corresponding to an SNR ≈ 80 dB, for "Barbara", the gain in ISNR for structured illumination is more than threefold, 7.12 in comparison to the deconvolution with 1.68 (for the same amount of data). In the case of "test pattern" too, the gain is significant, 6.1 for structured illumination and 3.5 for deconvolution. In both cases, there is an improvement in the SSIM index from 0.75 for the deconvolution to 0.85 for the structured illumination. The quantitative improvement of the classical approach in all cases is easily explained by the use of more than double the number of images than with our approach, reducing the influence of noise.

When the noise level is quite high (SNR ≈ 20 dB), the results may seem quite disappointing. Indeed, although an improvement over raw data is visible, the performances of SI are not much higher than those of the deconvolution. This can be explained by the fact that inside the support of the optical transfer function, the white noise corrupts the aliased high frequencies coming from the modulation. In other words, even if high frequencies are aliased inside the support of the transfer function, because their level is low, they are indistinguishable from the noise. As described in Sec. II-E, the proposed approach allows to compute the uncertainty about the estimated image through the standard deviation of the posterior law. Fig. 11 illustrates the estimated uncertainty with the "test pattern" with SNR ≈ 80 dB, while Tab. I reports the mean value. In this experiment, the uncertainty σ has an average value of ≈ 15 for a signal value between 0 and 255. A very striking fact is the structure of the uncertainty (Fig. 11a) in relation to the illumination pattern (Fig. 11b). Clearly, there is a strong correlation between the illumination and the uncertainty patterns: a higher illumination of the sample, i.e. more photons being shone on and emitted from the sample, is associated with a lesser uncertainty on the reconstruction, and conversely when there are less photons, the uncertainty is higher. While this is true when working with simulated data with a stationary Gaussian noise model, it may differ in the case of a non stationary noise model like Poisson.

As can be seen from Tab. I, the use of more than double images makes the classical SI approach slightly less sensitive to noise with a reduced uncertainty. 2) Hyper-parameters estimation: Fig. 12 illustrates the Markov chains and histograms of the noise and image power for the "Barbara" experiment with SNR ≈ 80 dB. The chains look similar for other experiments.

Numerical estimation results are reported in Tab. II. For the parameter γ n , the estimations are very close to the true values in all the cases, γ n = 0.99 and γ n = 0.0029 for true values equal to 1 and 0.03 respectively. The uncertainties, provided by the standard deviation of the sample, are small, indicating that the data are informative enough about this parameter. The uncertainty in the case of a SNR ≈ 20 dB is much lower than in the 80 dB case. This may be interpreted as the increase in noise power introducing more information about γ n value.

As to γ f , since the true value is not known, the error value

(MSE) e(γ f ) = f -f (γ f ) 2 / f 2 (39)
is computed for several values of γ f , ranging from 10 -4 to 10 in logarithmic scale, to determine the optimal one for image reconstruction. Results are reported in Tab. II. Globally, the estimated γf is in the same range as the optimal one. The uncertainty is small as in the γ n case. Consequently, the data are informative enough to automatically tune the trade-off between the likelihood and the spatial regularity. The Markov chains and histograms of the parameters are illustrated in Fig. 12 for "Barbara" with an SNR ≈ 80 dB. Figs. 12a and12d illustrate the full chains from the initialization on. For both parameters, the burn-in time is approximately 10 samples. Figs. 12b and12e illustrate the chains from sample 20 to the end. There is a good exploration of the state space around a stable value corresponding to the computed mean. The dispersion of the sample corresponds to the uncertainty. Figs. 12c and 12f are histograms of the sample after the burn-in time. The histograms are concentrated around the stable value and their width corresponds to the uncertainty. These are representations of the marginal posterior law for these parameters. (f) Fig. 12. Chains of hyper-parameters for "Barbara" with SNR ≈ 80 dB. The symbol × denotes the initialization. The x-axis of the chains is iteration and y-axis is the precision value (square value of image unit). The x-axis of the histograms is the precision value and y-axis is the count of samples.

3) Recovering still higher frequencies with eight images: As noted in section III, our approach allows the use and explotation of any modulation strategy. This not only allows to reduce the number of necessary images from 9 to 4, but also allows tiling the space with additional data to reconstruct additional high frequencies. This idea is illustrated in Fig. 13. To extend the first ring, we need only use 4 additional LR images to reconstruct the full ring as explained in Sec. III. The only requirement is the ability to illuminate the sample with a higher frequency modulation, with saturated excitation for example [START_REF] Heintzmann | Saturated patterned excitation microscopy with twodimensional excitation patterns[END_REF]. The required number of images is then decreased to 8 (4+4) instead of 21 with the classical SI approach.

The results of this experiment are simulated in Fig. 14.

f x f y Fig.

B. Real data

We have tested our method on real fluorescence microscopy data. Fig. 16 is an example of image reconstruction through a SIM device built in our laboratory. The optical setup is based on a commercial microscope body (Olympus IX81). A schematic diagram of the system is presented in Fig. 15. Fig. 16a is one example of LR fluorescent images obtained with a standard wide field microscope [START_REF] Rodríguez | Axial coding in full-field microscopy using three-dimensional structured illumination implemented with no moving parts[END_REF]. The excitation laser light wavelength is 488 nm. The excitation filter is a band pass 488 ± 10 nm (Chroma Z488/10x) and the emmision filter is a long pass 500 nm (Chroma HQ500LP). The dichroic is a T495lpxr. Imaged sample in Fig. 16 are actin filaments labeled with quantum dots (QDs) and Figs. 17 and 18 are QDs is solution spread over a cover slip.

The standard illumination is replaced by a fringe projection setup to produce the amplitude modulation. Laser light was scrambled by a rotating diffuser to reduce spacial coherence (to have a better optical cut and reduce the influence of the out of focus light). The fringe pattern is generated by coherent illumination of a spatial light modulator (SLM) displaying a periodic structure [START_REF] Rosen | Longitudinal spatial coherence applied for surface profilometry[END_REF]. We selected only the +1 and -1 diffraction orders by beam blocking the zero and superior orders. Both beams were focused in the back focal plane of the objective lens which made them interfere to form the illumination pattern with lateral sinusoidal structure and extending roughly 350 nm in the axial direction. Parameters of the illumination were estimated using cross-correlation techniques described in [START_REF] Gustafsson | Threedimensional resolution doubling in wide-field fluorescence microscopy by structured illumination[END_REF], [START_REF] Mandula | Patterned excitation microsopy[END_REF]. A 16-bits camera was used with a ×100 zoom and a pixel size of 8 × 8 µm as pixel size.

Fig. 16b is the result of our reconstruction algorithm when using 9 LR images, 3 LR images in each orientation as described in Sec. I-C. The reconstruction Fig. 16c is the result obtained when using only 4 LR images extracted from the data set.

The results clearly show an improvement both in resolution and contrast. Filaments are more visible than in standard microscopy and small details are now visible. Structured illumination clearly improves spatial resolution. The results obtained with 9 and 4 images are very similar, demonstrating the capability of the proposed approach to effectively decrease the required number of LR images. The differences between Figs. 16b and 16c are mostly due to contrast.

Figs. 17 and 18 demonstrate the effectiveness of our approach also when used on different types of samples, like these images of quantum dots (QDs) on a flat surface. Note that here we are not seeing individual QDs but rather aggregates. The differences between the classical SIM approach 17b and our SIM method 17c are minimal. On Fig. 18 and Fig. 19, the gain in resolution brought in by structured illumination can be appreciated by the fact that the QD aggregates that are not resolved by standard wide-field microscopy are well separated by SIM (Fig. 19c), irrespective of using 9 or just 4 images.

V. CONCLUSION

This paper has presented a novel approach for image reconstruction in Structured Illumination Microscopy. The proposed algorithm allows overcoming several drawbacks of existing approaches like automatic parameter estimation, reduction of necessary images and quantitative evaluation of the uncertainty. This is done through an inverse problem approach conducted within the Bayesian framework.

As illustrated in Sec. IV-B, the reconstruction is very sensitive to the modulation parameters used in the algorithm which may constitute a serious difficulty in practice. Several works already exist to address this point [START_REF] Schaefer | Structured illumination microscopy: artefact analysis and reduction utilizing a parameter optimization approach[END_REF]- [START_REF]Phase-shift estimation in sinusoidally illuminated images for lateral superresolution[END_REF]. One solution would be the joint estimation of these parameters together with the image and hyper-parameters, all within the same framework. This approach would lead to a myopic or semiblind image reconstruction [START_REF] Orieux | Bayesian estimation of regularization and PSF parameters for Wiener-Hunt deconvolution[END_REF].

Other perspectives are the use of different models for the noise, like Poisson, often encountered in microscopy. Different and more edge-preserving models for the image can also be considered like TV [START_REF] Chantas | Variational Bayesian image restoration with a product of spatially weighted total variation image priors[END_REF] or Huber-like [START_REF] Huber | Robust Statistics[END_REF] potentials. N (m f , Σ f ), where

m f = Σ -1 f A t Σ -1 n g (43)
Σ f = A t Σ -1 n A + R -1 f . ( 44
)
If g and m are drawn by g ∼ N (g, Σ n ) and m ∼ N (0, R f )

which means that data and prior mean are perturbed with their respective prior law, then the minimum is defined by

f = arg min f ( g -Af ) t Σ -1 n ( g -Af ) + (f -m) t R -1 f (f -m) (46) = Σ -1 f A t Σ -1 n g + R f m (47)
following a Gaussian law with mean m f and covariance Σ f . Proof: The variable f is Gaussian as it is a linear combination of two Gaussian random variables. The mean and covariance are calculated below. From (47), we have

E[ f ] = Σ -1 f E A t Σ -1 n g + R -1 f m = Σ -1 f A t Σ -1 n E[ g] + R -1 f E[ m] = Σ -1 f A t Σ -1 n g = m f . (48)
From ( 47) and (48), we have

f -E[ f ] = Σ f A t Σ -1 n ( g -g) + R -1 f m .
So, the covariance matrix is

V[ f ] = E f -E[ f ] f -E[ f ] t = Σ f A t Σ -1 n E ( g -g)( g -g) t Σ -1 n A+ R -1 f E m m t R -1 f Σ f = Σ f A t Σ -1 n E δ g δ t g Σ -1 n A+ R -1 f E δ f δ t f R -1 f Σ f = Σ f A t Σ -1 n Σ n Σ -1 n A + R -1 f R f R -1 f Σ f = Σ f A t Σ -1 n A + R -1 f Σ f = Σ f Σ -1 f Σ f = Σ f .
Then f follows the target law N (m f , Σ f ).

Fig. 1 .

 1 Fig. 1. Amplitude modulation of a signal. The original pattern (solid black) is centered around 0. The replications, due to modulation (light gray), are centered about the modulation frequencies fm and -fm. Outside the support of the transfer function, the frequencies are null, represented in dashed line. With the modulation, High-frequencies from fc to 2fm are inside the support through the modulation.

Fig. 3 .

 3 Fig. 3. Modulation of a 1D signal. Because of the conjugacy, only halffrequency planes are shown.

Fig. 4 .

 4 Fig. 4. Modulation of a 2D signal. (a) Only two unknowns (U 0 , U 1 ) are present in the surface f 2 x + f 2 y ≤ fc except in the gray parts where three unknowns are present; (b) The gray areas represent the frequencies where the likelihood alone is not sufficient to resolve the ambiguity in the final reconstruction.

Fig. 5 .

 5 Fig. 5. True images used for the study (in arbitrary units).

Fig. 6 .

 6 Fig.6. Comparison of reconstruction quality on "Barbara" with SNR ≈ 80 dB. The x-axis of slice images gives the pixel number. The improvement achieved by structured illumination over deconvolution is clearly noticeable. The result of our SI approach with 4 images is very close to the classical SI approach with 9 images.

Fig. 7 .

 7 Fig. 7. Spectrum of "Barbara". Figures are the logarithm of the absolute value of FFT-2D.

Fig. 8 .

 8 Fig. 8. Empirical estimation of circular power spectral density. The x-axis is the reduced frequency. The y-axis is the square value of the absolute value of Fourier transform in dB.

Fig. 10 .

 10 Fig. 10. Spectrum of "test pattern". Figures are the logarithm of the absolute value of FFT-2D.

Fig. 9 .

 9 Fig. 9. Comparison of reconstruction quality on "test pattern". The improvement achieved through Structured Illumination with 9 or 4 images is clear.

Fig. 11 .

 11 Fig. 11. Uncertainty about the estimate for simulated data and a stationary Gaussian noise model. Fig. 11a is the posterior standard deviation (in image unit), see Eqn. (32) of the "test pattern" with SNR ≈ 80 dB. Fig.11bis the mean of all illumination patterns (without unit). Clearly, there is less uncertainty at locations where there is more illumination. Both axis are pixel numbers.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Illustration of reconstruction improvement through additional high frequencies measurements. The images on the bottom row are the spectra of the images on the upper row (both axis in reduced frequency).

Fig. 19 .

 19 Fig. 19. Profile of a line in Fig. 18. (a) zoom on QD aggregates in the widefield image. (b) zoom on QD aggregates in the SIM reconstructed image with 9 LR images. (c) profile of the lines through the QD aggregates widefield image (WFM) and the SIM reconstruction with 9 and 4 images. The gain in resolution can be appreciated. Field of view is 960x960 nm.

TABLE I NUMERICAL

 I EVALUATIONS OF PERFORMANCE. ISNR AND SSIM ARE THE GAIN PER LR IMAGES USED. DATA CORRESPOND TO LR IMAGES OBTAINED WITHOUT STRUCTURED ILLUMINATION (SI). [22] IS THE DECONVOLUTION, "PROPOSED" IS OUR RECONSTRUCTION WITH ONLY 4 IMAGES AND [13] THE CLASSICAL SI APPROACH WITH 9 IMAGES. QUANTITATIVE PERFORMANCE IS SLIGHTLY BETTER WITH [13] BECAUSE OF THE USE OF MORE THAN DOUBLE THE NUMBER OF IMAGES FOR THE SAME QUALITATIVE RESULT.

TABLE II ESTIMATION

 II OF HYPER-PARAMETERS. THE VALUES γ f * ARE THE OPTIMAL VALUES WITH RESPECT TO THE MSE (EQ. 39) OR THE SSIM (EQ. 38).

  [START_REF] Mandula | Patterned excitation microsopy[END_REF]. for estimation of a second ring. The black pattern (first ring) is obtained with the illumination presented in Fig.2. The gray pattern corresponds to the second ring. The gray dashed pattern is the conjugate of the gray pattern. Figs. 14b and 14e are the zoom and logarithm of the absolute value of DFT-2D of the "test pattern", with only the first ring. Figs. 14c and 14f are the results with the second ring, with 8 LR images. The SNR ≈ 80 dB for both experiments. The gain is clearly visible, in direct and in Fourier domains. Numerical evaluation provides ISNR ≈ 8.65 dB and SSIM ≈ 0.67 for 8 LR images in comparison to ISNR ≈ 6.1 dB and SSIM ≈ 0.65 for 4 LR images.

Coherent systems are addressed in Sec. IV-B

The term resolution is used in the sense of spectral content here and not in the sense of sampling frequency.

This project has been financed by ANR, Institut Pasteur and ESPCI. FO was funded through an ANR fellowship.

However, these kinds of models, because of the unknown normalisation constant of the law, do not allow a straightforward estimation of hyper-parameters.

APPENDIX A SIMULATION OF CONDITIONAL POSTERIOR LAW OF THE

IMAGE

The conditional posterior law of the image is Gaussian

where the mean and covariance matrix are

The difficulties are 1) the high dimension law (the number of pixels in HR image) and 2) the correlation present in Σ f that embed H. The proposed solution is based on the same approach as in [START_REF] Rue | Fast sampling of Gaussian Markov random fields[END_REF], [START_REF] Lalanne | Stochastic artificial retinas: algorithm, optoelectronic circuits, and implementation[END_REF] but the resolution is done with an iterative optimization algorithm instead of an explicit inversion of the matrix.

Proposition A.1: Consider a non-stationary Gaussian law


Sandrine Maljean-Dubois

Matthieu Wemaere

MAKING THE EFFECTIVE LINKS WITHIN THE CLIMATE REGIME COMPLEX: THE 2015 AGREEMENT

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

POLICY OPTIONS

The concept of regime complex has the advantage of highlighting what links the different regimes together as much as what divides them. Hence it offers an interesting Seen from this perspective, the Paris Accord should be viewed as the bedrock of this regime complex; otherwise the Accord and the whole regime complex are likely to be ineffective. To this end, the Accord should aim to fulfil two different but complementary objectives:

-on the one hand, a catalysing role to create a dynamic and contribute to raising the level of ambition in the other regimes that form part of the regime complex:

-on the other hand, to play the leading role in order to orchestrate climate governance, to strengthen coherence and ensure complementarity, ensure work is covered correctly and avoid duplication of effort etc.

PROMOTING MUTUAL SUPPORTIVENESS

international trade law. 1. In order to effectively address climate change and its adverse effects, States shall formulate, elaborate and implement international law relating to climate change in a mutually supportive manner with other relevant international law. 2. States in cooperation with relevant international organizations shall ensure that consideration of climate mitigation and adaptation will be integrated into their law, policies and actions at all relevant levels, as laid out in Article 3. 3. According to Article 8, States shall cooperate with each other to implement the inter-relationship principle in all areas of international law, whenever necessary (…)"