Edouard Ollier
email: edouard.ollier@ens-lyon.fr

V Viallon

Regression modeling on stratified data with the lasso

HAL is

Introduction

We consider the estimation of regression models when the population under study is stratified, as is standard in epidemiology and clinical research. For instance, when studying relapse after a primary breast cancer, it is now common to analyze various histological subtypes at once (Voduc et al., 2010;[START_REF] Rosner | Breast cancer risk prediction with heterogeneous risk profiles according to breast cancer tumor markers[END_REF]. In order to accurately estimate the risk of relapse according to cancer subtype, risk factors that interact with cancer subtype need to be identified, and the corresponding interactions need to be precisely described. In pharmacokinetics, it is often of interest to describe how parameters related to absorption and clearance depend on dosage and type of adjuvant. In [START_REF] Ollier | A SAEM algorithm for fused lasso penalized non linear mixed effect models: Application to group comparison in pharmacokinetic[END_REF], for example, non-linear mixed effect models are estimated on strata defined according to the treatment dosage and the type of adjuvant. In Section 4, we study the association between the expression of the epidermal growth factor receptor gene, and those of 44 other transcription factors at eight time points of a differentiation process. We use a data set described in [START_REF] Kouno | Temporal dynamics and transcriptional control using single-cell gene expression analysis[END_REF] where expression of these factors has been measured by profiling, for each time point, 120 different single cells. In this application, our aim is to describe how the association between the epidermal growth factor receptor gene and the other 44 transcription factors varies over these eight strata.

In all these examples, the general objective is to study the relationship between a response variable y ∈ IR and a vector of p ≥ 1 predictors x = (x 1 , . . . , x p ) ∈ IR p over K ≥ 1 strata defined according to a categorical covariate Z of primary interest, such as cancer subtype. A key objective is to determine how Z modifies the effect of x on y, that is to identify and describe interactions between predictors x and the categorical effect modifier Z [START_REF] Gertheiss | Regularization and model selection with categorial effect modifiers[END_REF]. Identifying and taking account of these interactions is important even when the assessment of categorical effect 1 arXiv:1508.05476v2 [math.ST] 8 Nov 2016 modification is not the main objective. In particular, estimating one model per stratum or one model on all the strata pooled together generally leads to overfitting or underfitting, respectively.

For simplicity, we mostly focus on linear regression models. Denote by β * k = (β * k,1 , . . . , β * k,p ) ∈ IR p , for any k ∈ [K] = {1, . . . , K}, the parameter vector describing the association between y and x in the model corresponding to stratum Z = k. For any j ∈ [p], denote by d j ∈ [K] the number of distinct values in the set {β * 1,j , . . . , β * K,j }. Further consider the partition {K

(1) j , . . . , K

(d j ) j } of [K] such that β * k 1 ,j = β * k 2 ,j if and only if (β * k 1 ,j , β * k 2 ,j ) ∈ K (d) j for some d ∈ [d j ].
The full description of how Z modifies the effect of x j on y relies on the identification of this partition. The total number of partitions of [K] is the Kth Bell number B K , with, for instance, B 5 = 52. Standard statistical test procedures would require the comparison of (B K ) p models for the identification of the p partitions, so they are not well-suited to fully describe how Z modifies the effect of x on y.

Penalized approaches have been advocated in this context, which can be seen as a special case of multi-task learning [START_REF] Evgeniou | Regularized multi-task learning[END_REF]. In particular, a version of the adaptive generalized fused lasso (Tibshirani et al., 2005;[START_REF] Ollier | A SAEM algorithm for fused lasso penalized non linear mixed effect models: Application to group comparison in pharmacokinetic[END_REF] has been shown to enjoy an asymptotic oracle property if Kp does not grow with the sample size n [START_REF] Gertheiss | Regularization and model selection with categorial effect modifiers[END_REF][START_REF] Oelker | Regularization and model selection with categorical predictors and effect modifiers in generalized linear models[END_REF]. If Kp is fixed then this method identifies partitions {K (1) j , . . . , K (d j ) j } for all j ∈ [p] with probability tending to one as n → ∞, under mild assumptions. However, theoretical results in a non-asymptotic framework are still lacking for this method. Results obtained by [START_REF] Sharpnack | Sparsistency of the edge lasso over graphs[END_REF] for generalized fused lasso estimates in the Gaussian means setting are not straightforward to extend to our context but they suggest that this method might only be able to identify the partitions under very particular settings. See Section 2.6 for more details.

An alternative strategy, which is standard in epidemiology, consists in first selecting a reference stratum, say the first one, then coding the other strata by indicators 1I(Z = k), and finally including them into the regression model, along with their interactions with the predictors. This corresponds to considering the decompositions β * k = β * 1 + δ * k , for all k ∈ [K], with δ * 1 = 0 p , the null vector in IR p . Then, the lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] can be used to identify null components in vectors β * 1 and δ * k , for k = 1. Of course, this strategy can generally identify only one element of each partition {K

(1) j , . . . , K (d j ) j

}. However, a natural question is whether it is sparsistent, that is whether it can identify the sets

S * 1 = {j ∈ [p] : β * 1,j = 0} and T * 1 = {(k, j) ∈ [K] × [p] : β * k,j = β * 1,j
} with high probability. If so, it could partially describe the effect of Z on the association between y and x, while results of [START_REF] Sharpnack | Sparsistency of the edge lasso over graphs[END_REF] suggest that the method based on the fused lasso generally fails to fully describe it.

In this article, we establish that the sparsistency of this basic approach depends on the choice of the reference stratum. We then present a refined approach which bypasses this arbitrary choice. Under linear regression models, we show that our proposal is sparsistent under conditions similar to those ensuring the sparsistency of the approach based on an optimal and covariate-specific choice of the reference stratum. In addition, our proposal can be implemented with available packages under a variety of models, at approximately the same computational cost as that of the basic approach.

Methods

Notations and setting

Methods are first presented in the linear regression model for ease of notation. Extensions to generalized linear models [START_REF] Mccullagh | Generalized Linear Models[END_REF] are briefly presented in Section 2.4.

For any positive integer m ≥ 1, define [m] = {1, . . . , m}. Let 0 m and 1 m be the vectors of size m with components all equal to 0 and 1 respectively, and let I m be for the (m × m) identity matrix. For any vector x = (x 1 , . . . , x m ) T ∈ IR m , let supp(x) = {j ∈ [m] : x j = 0} denote its support. We further set x q = ( j∈[m] |x j | q ) 1/q for any real number q ∈ (0, ∞), x ∞ = max j |x j | and x 0 = |supp(x)|, where |E| is the cardinality of the set E. For any set E ⊆ [m], let x E denote the vector of IR |E| with components (x j ) j∈E . For any real matrix M , let M j be its jth column, and let M E be the sub-matrix made of columns (M j ) j∈E . Further denote the smallest singular value of M by Λ min (M ). Finally, 1I(•) is the indicator function.

Denote the number of levels of variable Z, that is the number of strata, by K ≥ 1. Let n k be the number of observations in stratum k ∈ [K], and denote the total number of observations by

n = k∈[K] n k . For k ∈ [K], further denote the response vector in stratum k by y (k) ∈ IR n k . Similarly, let X (k) be the (n k × p) design matrix in stratum k. For all k ∈ [K], we assume that y (k) = X (k) β * k + ε (k) , with noise vector ε (k) = (ε (k) 1 , . . . , ε (k) n k ) T ∈ IR n k . Vectors β * k ∈ IR p
include the Kp parameters to be estimated.

Basic approach

A basic approach consists in picking a reference stratum for any j, say j . Most often in practice, j is chosen so that it does not depend on j; here we consider the most general version. Introduce

= ( 1 , . . . , p ) ∈ [K] p , µ * = (β * 1 ,1 , . . . , β * p,p ) and δ * k = β * k -µ * .
The basic approach relies on the decomposition

β * k = µ * + δ * k , for all k ∈ [K]
, with δ * j ,j = 0 [START_REF] Gertheiss | Regularization and model selection with categorial effect modifiers[END_REF]. Following the lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], estimates of µ * and δ * k , k ∈ [K], can be defined as

argmin µ,δ 1 ,...,δ K δ j ,j =0 for all j∈[p] K k=1 y (k) -X (k) (µ + δ k ) 2 2 2n + λ 1 µ 1 + K k=1 λ 2,k δ k 1 , (1)
for appropriate non-negative λ 1 and λ 2,k . As will be seen in Sections 2.5 and 2.6, conditions ensuring the sparsistency of this approach depend on the arbitrary choice of the vector of reference strata . As a matter of fact, the underlying model dimension is µ * 0 + k = β * k -µ * 0 , which depends on . The lowest dimension, and then the best possible performance, is attained if µ * j is a mode of the collection of values (0, β * 1,j , . . . , β * K,j ), mode(0, β * 1,j , . . . , β * K,j ), for all j ∈ [p]. For any j ∈ [p], such an optimal and covariate-specific reference stratum will be denoted by * j below. Because * = ( * 1 , . . . , * p ) is generally unknown, the corresponding optimal version of the basic approach cannot be implemented in practice.

Our proposal

Our proposal aims at bypassing the arbitrary choice of the vector of reference strata , while mimicking the optimal version of the basic approach. We consider an overparametrization involving

(K + 1)p parameters, β * k = µ * + γ * k for any k ∈ [K].
It generalizes the decomposition used in the basic approach in the sense that no coefficient is constrained to be zero here. For nonnegative values of λ 1 and λ 2,k 's, our proposal returns

( µ, γ 1 . . . , γ K ) ∈ argmin µ,γ 1 ,...,γ K K k=1 y (k) -X (k) (µ + γ k ) 2 2 2n + λ 1 µ 1 + K k=1 λ 2,k γ k 1 . (2)
Working with a large enough value for λ 2,r , for some r ∈ [K], is equivalent to constraining γ r = 0 p , and then µ = β r , and reduces to the basic approach with = (r, . . . , r). More generally, setting τ = (τ 1 , . . . , τ K ) with τ k = λ 2,k /λ 1 , and defining the shrunk and τ -weighted version of the median of (b 1 , . .

. , b K ) as WSmedian(b 1 , . . . , b K ; τ ) = argmin b (|b| + k∈[K] τ k |b k -b|)
, it is easy to see that µ j ∈ WSmedian( β 1,j , . . . , β K,j ; τ ). In other words, for any particular value of the λ 2,k /λ 1 ratios, our approach encourages solutions ( β 1 = µ + γ 1 , . . . , β K = µ + γ K ) with a sparse vector µ and sparse vectors of differences γ k = β k -µ, and with the overall effect of the jth covariate µ j defined as WSmedian( β 1,j , . . . , β K,j ; τ ). Moreover, working with a large enough value for λ 1 is equivalent to constraining µ = 0 p , and our approach then reduces to K independent lasso's. In contrast, working with large enough λ 2,k values is equivalent to constraining β k = µ for all k, and our approach then reduces to one lasso run on all the strata pooled together.

Rewriting as a lasso on a transformation of the original data

Our proposal reduces to the lasso on a simple transformation of the original data, just as the basic approach does. Set Y = (y (1) T , . . . , y (K) T ) T the vector containing the n observations of the response variable. For any k ∈ [K], introduce P (k) = {j ∈ [p] : k = j }, with j still denoting the reference stratum chosen for covariate j in the basic approach. Note that k |P

(k) | = (K -1)p and set X(k) = X (k) P (k) . Now introduce X =    X (1) X(1) /τ 1 . . . 0 . . . . . . . . . . . . X (K) 0 . . . X(K) /τ K    , X 0 =    X (1) X (1) /τ 1 . . . 0 . . . . . . . . . . . . X (K) 0 . . . X (K) /τ K    ,
Criteria to be minimized in (1) and (2) reduce to the lasso

1 2n Y -X θ 2 2 + λ 1 θ 1 , (3)
with X set to either X , for the basic approach, or X 0 , for our proposal, and θ a vector of IR Kp or IR (K+1)p as appropriate. Therefore, our proposal comes at almost no additional computational cost compared to the basic approach. This rewriting as a lasso extends to generalized linear models, Cox models, etc., and makes our proposal directly implementable using the glmnet package of [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF], for instance. Considering logistic models, set L logistic (y, z) = i∈[n] y i z i -log(1 + e z i ) for any y ∈ {0, 1} n and z ∈ IR n . The criterion of our proposal writes as the logistic lasso,-L logistic (Y , X 0 θ) + λ 1 θ 1 . In addition, the glmnet package can benefit from the sparse structure of X 0 .

Sparsistency

From now on we assume that β * j = mode(0, β * 1,j , . . . , β * K,j ) is uniquely defined for all j ∈ [p], for ease of notation. For any vector of reference strata = ( 1 , . . . , p ) ∈

[K] p , introduce the sets S = {j ∈ [p] : β * j ,j = 0} and T = {(k, j) ∈ [K] × [p] : β * k,j = β * j ,j }. Further define θ * = (µ * T , τ 2 γ * T ,2 , . . . , τ K γ * T ,K ) T ∈ IR Kp , with µ * ,j = β * ,j and γ * ,k = (β * k -µ * ) P ( ) k
. The basic approach is sparsistent if it identifies these two sets, or, equivalently, the set J = supp(θ * ), with high probability. Now, consider an optimal vector of reference strata * = ( * 1 , . . . , * p ) such that β * * j ,j = β * j for all j ∈

[p]. Define θ * 0 = (µ * T * , τ 1 γ * T 0,1 , . . . , τ K γ * T 0,K ) T ∈ IR (K+1)p where γ * 0,k = β * k -µ * * .
We will say our proposal is sparsistent if it identifies S * and T * , or, equivalently, the set J 0 = supp(θ * 0 ), with high probability. For the lasso to be sparsistent, a sufficient and almost necessary condition on the design matrix is the irrepresentability condition (Zhao & Yu, 2006;Wainwright, 2009). Consider the general formulation (3) of the lasso and denote by θ * the true value of the parameter vector to be estimated. Defining J * = supp(θ * ), the matrix X fulfills the irrepresentability condition, with respect to J * , if Λ min (X T J * X J * ) ≥ C min for some C min > 0 and max j / ∈J * (X T J * X J * ) -1 X T J * X j 1 < 1. Using the vector of reference strata in the basic approach, the irrepresentability condition of matrix X writes as (IC) , while in our approach, the irrepresentability condition of matrix X 0 writes as (IC) 0 :

(IC) Λ min (X T J X J ) ≥ C > 0 and c = max j / ∈J (X T J X J ) -1 X T J X j 1 < 1, (IC) 0 Λ min (X T 0J 0 X 0J 0 ) ≥ C 0 > 0 and c 0 = max j / ∈J 0 (X T 0J 0 X 0J 0 ) -1 X T 0J 0 X 0j 1 < 1.
The comparison of (IC) * and (IC) 0 is of particular interest. First, because θ *

* J * = θ * 0J 0 , we have X * J * = X 0J 0 and Λ min (X T * J * X * J * ) = Λ min (X T 0J 0 X 0J 0 ).
Second, the maxima in the definitions of c 0 and c * are taken over J c 0 and J c * , respectively, with

|J c 0 | = p + |J c * |. Indeed, columns corre- sponding to γ * 0, * j ,j , for j ∈ [p]
, are present in matrix X 0 while they are absent from X * . Moreover, the corresponding indexes belong to J c 0 since γ * 0, * j ,j = 0. Therefore, the only difference between (IC) 0 and (IC) * comes from the fact that c 0 ≥ c * , and (IC) 0 is only slightly stronger than (IC) * . Focusing on the case of balanced strata and orthogonal designs in each stratum, Lemma 2 in Section 2.6 states that these two conditions are identical in this particular case. It further explicitly relates them to the ratios τ k = λ 2,k /λ 1 and to the maximum level of heterogeneity that is allowed among the collections of values (β * 1,j . . . , β * K,j ), for all j ∈ [p]. In addition, they are shown to be generally weaker than (IC) for the choice = (r, . . . , r), for some r ∈ [K].

We can now state Theorem 1, according to which our proposal identifies S * and T * under nearly the same assumptions as those required by the optimal version of the basic approach if the * j 's were given in advance. More precisely, besides the fact that (IC) 0 is generally a little stronger that (IC) * , the only difference lies in the terms (λ

(1) 1 , β (1) min ) and (λ (0) 1 , β (0) min ), where K + 1 replaces K. Our result is a consequence of Theorem 1 in Wainwright (2009).
Theorem 1 Assume that the noise variables (ε

(k) i ) i∈[n k ],k∈[K]
are independent and identically distributed centered sub-Gaussian variables with parameter σ > 0. Further assume that n

-1/2 k X (k) j 2 ≤ 1 for all (k, j) ∈ [K] × [p]. Introduce τ 0 > 0 and set τ k = τ 0 (n k /n) 1/2 for all k ∈ [K]. If (IC) * holds then set γ 1 = 1 -c * . Further set γ 0 = 1 -c 0 if (IC) 0 holds. For η ∈ {0, 1}, introduce λ (η) 1 > 2 γ η min(1, τ 0 ) 2σ 2 log((K + η)p)) n 1/2 , β (η) min = λ (η) 1 (|S * | + |T * |) 1/2 C * + 4 σ C 1/2 * . If (IC) * holds then solutions θ * of (3) with X = X * and λ 1 = λ (0)
1 as above are such that, with probability at least 1 -4 exp(-a 0 nλ 1 ) for some constant a 0 > 0: (i) θ * is uniquely defined,

(ii) θ * J c * = 0 |J c * | , and (iii) θ * J * -θ * * J * ∞ ≤ β (0) min . If, in addition, | β * j | > β (0) min for all j ∈ S * and |β * k,j -β * j | > β (0)
min /τ k for all (k, j) ∈ T * , then J * , hence both S * and T * , are perfectly identified with probability at least 1 -4 exp(-a 0 nλ 2 1 ). If (IC) 0 holds then solutions θ 0 of (3) with X = X 0 and λ 1 = λ

(1) 1 as above are such that, with probability at least 1 -4 exp(-a 1 nλ 1 ) for some constant a 1 > 0: (i) θ 0 is uniquely defined, (ii)

θ 0J c 0 = 0 |J c 0 | , and (iii) θ 0J 0 -θ * * J * ∞ ≤ β (1) min . If, in addition, | β * j | > β (1)
min for all j ∈ S * and

|β * k,j -β * j | > β (1)
min /τ k for all (k, j) ∈ T * , then J 0 , hence both S * and T * , are perfectly identified with probability at least 1 -4 exp(-a 1 nλ 2 1 ).

This result especially confirms that it is harder to identify T * than S * , in the sense that heterogeneities have to be at least

|β * k,j -β * j | > (n/n k ) 1/2 β (η) min /τ 0 for (k, j) ∈ T * , while | β * j | has only to be greater than β (η) min for j ∈ S * .

The orthogonal and balanced case

It is instructive to inspect in more detail the simple setting where

n k = n/K and (X (k) T X (k) )/n k = I n k for all k ∈ [K]
. This orthogonality assumption does not make matrices X and X 0 orthogonal and is therefore not sufficient for the irrepresentability condition.

For any vector of reference strata and all j ∈

[p], define K * ,j = {k ∈ [K] : β * k,j = β * j ,j }. Set D ,0 = max j / ∈S |{k ∈ [K] : β * k,j = β * j ,j }| if S = [K] and 0 otherwise and D ,1 = max j∈S |{k ∈ [K] : β * k,j = β * j ,j }| if S = ∅ and -∞ otherwise. For all k ∈ [K], we set τ k = τ 0 K -1/2 , for some τ 0 > 0.
Lemma 2 The matrix X fulfills the irrepresentability condition if and only if

(sIC) 0 ≤ K 1/2 K -2D ,1 < τ 0 < K 1/2 D ,0 •
The matrix X 0 fulfills the irrepresentability condition if and only if

(sIC) 0 0 ≤ K 1/2 K -2D * ,1 < τ 0 < K 1/2 D * ,0 •
Conditions (sIC) 0 and (sIC) * in Lemma 2 are identical. In the orthogonal and balanced case, the two sets of assumptions required by our proposal and the optimal version of the basic approach to identify the sets S * and T * are therefore identical, except for the terms (λ

(1) 1 , β (1) min ) and (λ (0) 1 , β (0)
min ) where K + 1 replaces K, as in Theorem 1 above. In addition, (sIC) * , or equivalently (sIC) 0 , imposes that 2D * ,1 + D * ,0 < K. In particular, if D * ,1 = D * ,0 = D * , this implies that D * < K/3. The irrepresentability conditions (sIC) * and (sIC) 0 are therefore directly related to the maximum level of heterogeneity among the values (β * 1,j , . . . , β * K,j ). Similarly, (IC) imposes 2D ,1 + D ,0 < K, which is generally a stronger constraint. For simplicity, consider the situation where = (r, . . . , r) for some r ∈ [K], which is a common choice in practice. Without loss of generality, set r = 1. Then (IC) entails that, for each j ∈ [p], the effect of the jth covariate on most strata is β * 1,j , while (sIC) * and (sIC) 0 only entail that, for each j ∈ [p], the effect of the jth covariate on most strata is mode(0, β * 1,j , . . . , β * K,j ).

Moreover, if = (1, . . . , 1) and {1} / ∈ ∩ j∈[p] K * * ,j then we have T = T * and, possibly, S = S * : the identification of T and S is both less interesting and less likely and our proposal should be preferred over the basic approach.

To recap, our non-asymptotic analysis shows that the partial description of the categorical effect modification due to Z through decompositions of the type β * k = µ * + γ * k is guaranteed only when the level of heterogeneity is not too high, that is when the number of nonzero components in vectors (γ * 1,j , . . . , γ * K,j ), for j ∈ [p], is not too high. In particular, the lowest level of heterogeneity is attained for the choice µ * j = mode(0, β * 1,j , . . . , β * K,j ). Our proposal is able to target this optimal decomposition and is sparsistent under nearly the same assumptions as those required for the optimal version of the basic approach.

Connection with the generalized fused lasso

The fact that the level of heterogeneity must be low to ensure the sparsistency of our proposal as well as the sparsistency of the optimal version of the basic approach has connections with other results in the literature. Consider for instance the approach mentioned in Section 1, which was proposed by [START_REF] Gertheiss | Regularization and model selection with categorial effect modifiers[END_REF]; see also [START_REF] Oelker | Regularization and model selection with categorical predictors and effect modifiers in generalized linear models[END_REF] and [START_REF] Ollier | A SAEM algorithm for fused lasso penalized non linear mixed effect models: Application to group comparison in pharmacokinetic[END_REF]. This is based on a fusion penalty that encourages similarities among solutions

β k ∈ IR p , k ∈ [K].
More precisely, for appropriate non-negative λ 1 and λ 2 's, it returns estimators β k defined as minimizers of the criterion

k∈[K] y (k) -X (k) β k 2 2 2n + λ 1 k∈[K] β k 1 + λ 2 k 1 <k 2 β k 1 -β k 2 1 . (4)
This criterion is that of the generalized fused lasso, where the graph used in the penalty is made of p cliques of size K [START_REF] Ollier | A SAEM algorithm for fused lasso penalized non linear mixed effect models: Application to group comparison in pharmacokinetic[END_REF]. The jth clique corresponds to the jth predictor and connects all the components in (β 1,j , . . . , β K,j ): the K(K -1)/2 differences |β k 1 ,j -β k 2 ,j |, for k 1 < k 2 , appear in the penalty term.

Both the basic approach and our proposal are related to generalized fused lasso estimates too. In particular, consider criterion (1) with the optimal choice * for the vector of reference strata. It can be seen as a version of a generalized fused lasso, with a graph made of p star-graphs of size K instead of p cliques: for each j ∈ [p], only the K -1 differences |β k,j -β * j ,j | appear in the penalty term, for k = * j and * j ∈ [K] fixed. Non-asymptotic analyses of the sparsistency of generalized fused lasso estimates are scarce in the literature. [START_REF] Sharpnack | Sparsistency of the edge lasso over graphs[END_REF] study them in the normal means setting, which can be seen as a special case of the stratified linear regression considered here. [START_REF] Sharpnack | Sparsistency of the edge lasso over graphs[END_REF] establish that generalized fused lasso estimates are sparsistent only if the graph used in the fused penalty is in good agreement with the true structure of the vector of parameters; see also [START_REF] Qian | On stepwise pattern recovery of the fused Lasso[END_REF]. Although it is not straightforward to extend to our case, these results suggest that estimates derived from (4) can only be sparsistent if the level of heterogeneity is not too high.

Our results precisely quantify the maximum level of heterogeneity above which a version of generalized fused lasso estimates, based on star-graphs, can attain sparsistency in stratified regression, in the balanced and orthogonal case. Because star-graphs are less connected than cliques, it is likely that sparsistency for clique-based estimates, such as those minimizing criterion (4), requires an even lower maximum level of heterogeneity. That being said, sparsistency for clique-based estimates refers to the full identification of the partitions {K

(1) j , . . . , K (d j ) j }, for j ∈ [p].
For the basic approach, its optimal version and our proposal, sparsistency refers to the identification of one element of this partition only, say K

(1) j , and its complementary [K] \ K

(1) j . In the asymptotic regime, assuming that Kp is fixed and n k /n → ρ k for some ρ k ∈ (0, 1), for all k ∈ [K], oracle properties have been derived for adaptive versions of clique-based estimates under mild assumptions [START_REF] Gertheiss | Regularization and model selection with categorial effect modifiers[END_REF]: in particular, no assumption regarding the level of heterogeneity is required to ensure perfect recovery of the full partition {K

(1) j , . . . , K (d j ) j }, for all j ∈ [p]. Similar results are easily derived for adaptive versions of the basic approach for instance. In view of (3), we can apply Theorem 2 of Zou (2006) to show that an adaptive version of the basic approach enjoys an oracle property too, without having to assume any irrepresentability condition. Here again, no assumption regarding the maximum level of heterogeneity is required, but the identification of only one element of the partition is guaranteed for all j.

To recap, clique-based estimates are optimal and should be preferred over our proposal or the basic approach in the asymptotic regime, assuming that Kp is fixed and n k /n → ρ k for some ρ k ∈ (0, 1), for all k ∈ [K]. In a non-asymptotic setting, results for clique-based estimates are still lacking, while conditions ensuring the sparsistency of our proposal and the basic approach are established in the present article. In the following simulation study, we especially compare the clique-based strategy and our proposal on finite samples.

Simulation Study

Theorem 1 states that our proposal and the optimal version of the basic approach perform similarly with regard to the identification of S * and T * for appropriate values of λ 1 and τ 0 , under technical assumptions on the design matrices. The main objective of this simulation study is to assess the empirical performance of our proposal under general designs, and for λ 1 and τ 0 selected by 5-fold cross-validation. Comparisons are made with the basic approach with the choice = (1, . . . , 1) as well as an optimal choice * . Clique-based estimates are also considered.

We set K = 20 and take n k ∈ {10, 50, 100} and p ∈ {20, 100, 500}. For each k ∈ [K], rows of the design matrix X (k) are drawn from an N (0 p , Σ) distribution, with Σ the (p × p) Toeplitz matrix with element (i, j) equal to 0.5 |i-j| . We then randomly select a subset P 0 ⊂ [p] of size 20 and set β * k,j = 0 for all j / ∈ P 0 and k ∈ [K]. As for the values (β * 1,j , . . . , β * K,j ) for j ∈ P 0 , we consider four levels of heterogeneity d H . More precisely, for any given d H ∈ {1, 3, 6, 9}, we set β * k,j = 1 for k > d H , and β * k,j = 1 + δ * k,j for k ≤ d H for 10 indexes j randomly selected in P 0 . For the other 10 indexes in P 0 , we set β * k,j = 1 for k ≤ d H , and β * k,j = 1 + δ * k,j for k > d H . We further consider two cases for the δ * k,j values: they are either constantly set to K 1/2 or drawn from the uniform distribution on [K 1/2 /2, 2K 1/2 ] and then multiplied by ±1 (with probability 1/2). When the δ * k,j 's are constant, the collection of values (β * 1,j , . . . , β * K,j ) for each covariate j ∈ P 0 is made of two groups of distinct values, of sizes K -d H and d H . This situation should favor clique-based estimates. When δ * k,j is random, the collection of values (β * 1,j , . . . , β * K,j ) for each covariate j ∈ P 0 is made of d H +1 groups, one of size K -d H , and the other d H of size 1. This situation should favor our proposal and the optimal version of the basic approach. Observations of the response variable are then generated according to y

(k) = X (k) β * k + ε (k) , with each component of ε (k) drawn from an N (0, σ 2 ) distribution. The variance σ 2 is set to k∈[K] X (k) β * k replicates of data (X (k) , y (k) ), k ∈ [K].
Our results correspond to averages over these 50 replicates; see Figure 1.

In all the configurations, β * 20,j = mode(0, β * 1 , . . . , β * 20,j ) for all j ∈ [p]. We then set * j = 20 for all j ∈ [p] for the optimal version of the basic approach. On the other hand, β * 1,j = mode(0, β * 1 , . . . , β * 20,j ) for all j ∈ P 0 . Under this setting, which is of course extreme, the comparison between the results obtained using either or * for the reference strata allows a precise description of the impact of the reference stratum on the performance of the basic approach. Top panels of Figure 1 present results regarding the identification of the sets

T * P 0 = {(k, j) ∈ [K] × P 0 : β * k,j = β * * j ,j }
for the optimal version of the basic approach, our proposal and clique-based estimates, and that of the set

T * 1,P 0 = {(k, j) ∈ [K] × P 0 : β * k,j = β * 1,j }
for the basic approach. Here, we only consider covariates in P 0 because they are those that are the most differently accounted for by the four approaches we compare.

In the constant δ * k,j case, our empirical results clearly illustrate our theoretical ones. First, our proposal performs nearly as well as the optimal version of the basic approach. Second, the lower d H , the better they perform, as expected since d H = D 1 and D 0 = 0 here. Third, it is more difficult to recover T * 1,P 0 than T * P 0 , which is also expected since D

(1) 1

= K -d H > D 1 .
For the same reason however, as d H increases, T * 1,P 0 is easier to recover. A nice symmetry appears between the performance of our proposal and the optimal version of the basic approach on the one hand and the basic approach using = (1, . . . , 1), on the other hand. In addition, the recovery of the sets T * P 0 and T * 1,P 0 is only marginally affected by p. Results in the random δ * k,j case are mostly consistent with those in the constant case, except for the recovery of T * 1,P 0 which is mostly due to the fact that T * 1,P 0 = [K -1] × P 0 irrespective of d H in this random case. Finally, the clique-based strategy performs similarly to, or a little worse than, our proposal with respect to this criterion.

The bottom panels of Figure 1 present results regarding log( et al., 2017), which is a measure of the prediction error. Overall, our proposal and the optimal version of the basic approach perform similarly. They both outperform the basic approach, especially in the random δ * k,j case. The clique-based strategy outperforms our proposal in the constant δ * k,j case if d H is high enough, but only when the n k /p ratio is not too small. In the random δ * k,j case, our proposal clearly outperforms the clique-based strategy when p = 500, and more generally as d H increases and/or the n k /p ratio decreases. These results suggest that the clique-based strategy might not be able to fully account for, nor benefit from, the true structure in a high-dimensional setting. They also suggest that our proposal is better suited for this high-dimensional setting, as long as heterogeneity is not too high.

k∈[K] X (k) (β * k -β k ) 2 2 /n) (Dalalyan

Application on single-cell data

We analyse data describing mielocytic leukemia cells undergoing differentiation to macrophage. Expression levels of 45 transcription factors are measured at K = 8 distinct time points of this differentiation process (H0, H1, H6, H12, H24, H48, H72 and H96). Each time point defines a stratum where data on n k = 120 single cells are available. This data set is described in [START_REF] Kouno | Temporal dynamics and transcriptional control using single-cell gene expression analysis[END_REF]. In this application, the main objective is to determine how associations among the 45 transcription factors vary over time. [START_REF] Kouno | Temporal dynamics and transcriptional control using single-cell gene expression analysis[END_REF] focus on marginal associations and use univariate analyses while graphical models, which describe conditional associations, might be better suited. Their inference can be reduced to the identification and description of the neighborhood of each covariate (Meinshausen Bühlmann, 2006). Here, as a first step, we study how the neighbor-

* P 0 = {(k, j) ∈ [K] × P 0 : β * k,j = β * * j ,j }
for the other approaches. The higher, the better. The bottom three panels illustrate the prediction error; the lower, the better. Results are presented for both the constant and random δ * k,j cases. All results correspond to averages over 50 replicates in each configuration. Solid line: our proposal. Dotted line: optimal version of the basic approach. Dash-dot line: basic approach. Dashed line: clique-based approach. hood of one particular transcription factor, EGR2, varies over time. Towards this end, we consider stratified linear regression models that relates EGR2 to the other p = 44 factors on the K = 8 strata. Expression levels of EGR2 are centered within each stratum, and no intercept is included in the models. Then, parameters of interest are vectors β * 1 , . . . , β * 8 , where β * k ∈ IR 44 describes the association between EGR2 and the p = 44 transcription factors at the kth time point. We compare estimates returned by our proposal and four competitors. The basic approach is considered with two distinct choices for the reference stratum: we set it to either H0 or H96 for each covariate. We further consider the clique-based strategy. Finally, given the ordinal nature of the strata in this particular example, the variant based on chain graphs [START_REF] Gertheiss | Regularization and model selection with categorial effect modifiers[END_REF] can be seen as the reference method. We include it as well, even if our main objective in this illustrative application is to compare the other four approaches, which do not account for this additional information. For each approach, regularization parameters are selected by 5-fold cross-validation. Results are presented in Figure 2. For each method, the estimates correspond to a 44 × 8 matrix of the form ( β 1 , . . . , β K ) with K = 8 and β k ∈ IR 44 for any k ∈ [K]. Therefore, for any j ∈ [p], the jth row of each matrix corresponds to estimates ( β 1,j , . . . , β K,j ) of the effect of the jth covariate over the K time points, returned by the corresponding approach. This heat map representation allows an easy comparison of the pattern identified among the effects ( β 1,j , . . . , β K,j ) of each covariate across the 8 strata. Our first objective is to illustrate the impact of the reference stratum when using the basic approach. Considering for instance the association between EGR2 and MYB, most approach identify the same pattern: the association is constant between H1 and H96, while it is lower, or even null, at H0. However, setting the reference stratum to H0, the basic approach does not detect any heterogeneity. These results are consistent with what is expected if the true pattern is the one identified by the other approaches: the basic approach used with H0 as the reference stratum is unlikely to identify the heterogeneity if it occurs at H0. Now consider the association with ELK1. Our proposal, the basic approach with H0 as the reference stratum and the strategy based on chain graphs all suggest the absence of association between EGR2 and ELK1 on the time interval H0 and H72, and a positive association at H96. If the reference stratum is set to H96, the basic approach suggests a quite different pattern, which is again expected if the true pattern is the one returned by the other approaches. These results confirm that the reference stratum is critical for the basic approach. They further suggest that our proposal is able to identify appropriate covariate-specific reference strata.

The comparison of the patterns returned by our approach and the two fused lasso strategies mostly highlights that the clique-based strategy identifies fewer heterogeneities and that strategy based on chain graphs returns smoother patterns. Again, these results were expected given the connectivity of the clique and chain graph, respectively. Prediction error was evaluated by double 5-fold cross-validation. Among the approaches that do not account for the ordering of the strata, the best prediction error is obtained with our proposal, while the worst is 1.8% higher and is obtained with the clique-based strategy. The chain graph strategy leads to an improvement of 1.8% compared to our approach.

Two main conclusions can be drawn. When data come from several strata of the population and no information is available regarding which strata are likely to share similar effects, our proposal is a competitive approach. When additional information is available, as in this particular application where strata are naturally ordered, accounting for it can be beneficial.

Discussion

After submitting a first version of this work, we became aware of concurrent work by [START_REF] Gross | Data Shared Lasso: A novel tool to discover uplift[END_REF] where the authors introduce similar ideas. They apply it, in particular for the uplift problem in clinical research where the objective is to find sub-populations in a randomized trial for which an intervention is beneficial. In addition, there has been a recent line of works on penalized approaches aiming at identifying interactions, not necessarily between a categorical covariate and other predictors [START_REF] Lim | Learning interactions via hierarchical group-lasso regularization[END_REF][START_REF] Radchenko | Variable selection using adaptive nonlinear interaction structures in high dimensions[END_REF]. In this general context, strong hierarchy is often imposed: whenever an interaction between two variables is included in the model, the corresponding main effects are included too. However, this strong hierarchy is not desirable in our setting, where a coefficient can be nonzero in only one of the strata [START_REF] Gross | Data Shared Lasso: A novel tool to discover uplift[END_REF]. Moreover, when applied in our setting, these approaches can be seen as versions of the basic approach (1) based on extensions of the L 1 -norm penalty. In particular, a reference stratum has to be chosen as a first step, and the performance of these approaches would depend on this arbitrary choice. For instance, [START_REF] Radchenko | Variable selection using adaptive nonlinear interaction structures in high dimensions[END_REF] establish assumptions under which their approach is sparsistent. When applied with the reference strata , their main condition on the design matrix is very similar to (IC) , which is generally stronger than (IC) * and (IC) 0 . Then, these approaches could benefit from the ideas we developed in this article.

Our proposal is based on an overparametrization, which naturally raises the question of identifiability. We refer to [START_REF] Gross | Data Shared Lasso: A novel tool to discover uplift[END_REF] for some discussion. We shall add that there is no identifiability issue under the conditions of Theorem 1. If these conditions do not hold, and in particular if (IC) * is not fulfilled, even the optimal version of the basic approach is not sparsistent, and the identifiability issue related to the overparametrization is secondary.

Prediction bounds for our proposal can be derived under the conditions presented in this work. But weaker conditions might be sufficient, following recent work studying the lasso for correlated designs [START_REF] Dalalyan | On the prediction performance of the lasso[END_REF]. Another extension might concern the derivation of valid p-values or confidence intervals for the nonzero parameters identified by our proposal. Given its connection with the lasso, this post-selection inference might be derived by extending recent strategies proposed for lasso estimates [START_REF] Lee | Exact post-selection inference, with application to the lasso[END_REF].

We also plan to extend our proposal to other regression models, which is straightforward for a variety of models given its connection with the lasso. In particular, our proposal could easily be extended to stratified Cox models used in survival analysis when competing risks arise [START_REF] Rosner | Breast cancer risk prediction with heterogeneous risk profiles according to breast cancer tumor markers[END_REF], or to the conditional logistic models used in case-controls studies [START_REF] Reid | Regularization paths for conditional logistic regression: The clogitl1 package[END_REF]. The extension of clique-based estimates to other models is generally more computationally burdensome, partly because there is no proximal operator for the fused penalty. TIBSHIRANI, R. J., SAUNDERS, M., ROSSET, S., ZHU, J. & KNIGHT, K. (2005) Appendix: Supplementary Files 5.1 Technical details

PROOF OF LEMMA 1

Lemma 1 is established for matrix X 0 ; the proof for X follows from similar arguments and is omitted. Fix τ 0 > 0 and set

τ k = τ 0 K -1/2 for all k ∈ [K]. Recall that θ * 0 = (µ * * T , τ 1 γ * 0,1 T , . . . , τ K γ * 0,K T ) T , with µ * * j = β * * j ,j for j ∈ [p], and J 0 = {j ∈ [(K + 1)p] : θ * 0j = 0}.
For the sake of brevity, the proof is only presented in the case where S * = ∅ and T * = ∅, where

S * = {j ∈ [p] : µ * * j = 0} and T * = {(k, j) : β * k,j = µ * * j }. Setting T * k = {j ∈ [p] : (k, j) ∈ T * } and Σ (k) S * ,T k = X (k) T S * X (k) T k
we have

X 0J 0 =       X (1) S * X (1) T * 1 τ 1 . . . 0 . . . . . . . . . . . . X (K) S * 0 . . . X (K) T * k τ K       , (X T 0J 0 X 0J 0 ) =                   nI |S * | Σ (1) S * ,T * 1 τ 1 . . . . . . . . . Σ (K) S * ,T * K τ K Σ (1) T S * ,T * 1 τ 1 nI |T * 1 | τ 2 0 0 . . . . . . 0 Σ (2) T S *
Σ (K) T S * ,T * K τ K 0 . . . . . . 0 nI |T * K | τ 2 0                   .
(t, s) of C k is (C k ) (t,s) = -τ k /N * S * s if S * s = T * k,t and 0 otherwise (s ∈ [|S * |] and t ∈ [|T * k |]).
For k, ∈ [K] and k = , further introduce B k 1 ,k 2 , the matrix of size

|T * k 1 | × |T * k 2 | with element (t 1 , t 2 ) equal to τ k 1 τ k 2 /N * S * s if T * k 1 ,t 1 = T * k 2 ,
(X T 0J 0 X 0J 0 ) -1 = D C T C B .
Now, for any j / ∈ J 0 , the jth column X 0j of X 0 is of the form either (A) or (B):

(A) X 0j = (0 T , . . . , 0 T , X

(k 0 ) T j 0 , 0 T , . . . , 0 T ) T for some k 0 ∈ [K] and some j 0 / ∈ T * k 0 , (B) X 0j = (X (1) T j 0 , . . . , X (K) T j 0 ) T for some j 0 / ∈ S * . If X 0j is of form (A), then X 0j T X 0J 0 = ( dT , 0 T |T * 1 | , . . . , 0 T |T * K | ) ∈ IR |J 0 | with d ∈ IR |S * | and ds = n k 0 /τ k 0 if j 0 = S * s , 0 otherwise.
Therefore, if X 0j is of form (A), we have

X 0j T X 0J 0 (X 0J 0 T X 0J 0 ) -1 1 ≤ max j∈S * max k∈K * j n k τ k N * j (1 + / ∈K * j τ ). If X 0j is of form (B), then X T 0j X 0J 0 = (0 |S * | , d 1 , . . . , d K ) T ∈ IR |J 0 | , with d k ∈ IR |T * k | and d k,t = n k /τ k if j 0 = T * k,t , 0 otherwise.
In this case, we have

X 0j T X 0J 0 (X 0J 0 T X 0J 0 ) -1 1 ≤ max j / ∈S * / ∈K * j τ .
Moreover, under the setting considered in Lemma 1 we have n k = n/K for all k ∈ [K] and

τ k = τ 0 K -1/2 . Therefore, max j / ∈J 0 X 0j T X 0J 0 (X 0J 0 T X 0J 0 ) -1
1 < 1 if and only if assumption (sIC) 0 holds, which completes the proof of Lemma 1 for matrix X 0 .

VERSION OF THEOREM 1 IN THE SETTING OF ORTHOGONAL DESIGNS AND BALANCED

STRATA

Theorem 3 below is the version of Theorem 1 in the setting considered in Lemma 1, where n k = n/K and (X (k) T X (k) )/n k = I n k for all k ∈ [K]. We set D 0 = D * ,0 and D 1 = D * ,1 ; see Section 2.6 in the main text for the corresponding definitions. solutions of the characteristic polynomial of the matrix (X T 0J 0 X 0J 0 /n), p(λ) = det(X T 0J 0 X 0J 0 /n -λI |J 0 | ). Using the block structure of matrix (X T 0J 0 X 0J 0 /n -λI |J 0 | ), we get the following expression

p(λ) = (1 -λ) r 1 (τ -2 -λ) r 2 j∈S * ∩{∪ k T * k } λ 2 -λ(τ -2 + 1) + τ -2 1 - n -N * j n , with r 1 = |S * \{∪ k T * k }| and r 2 = |J 0 |-|S * |-|S * ∩{∪ k T * k }|. It follows that Λ min ((X T 0J 0 X 0J 0 /n)) ≥ min 1, τ -2 , 1 2 (τ -2 + 1) -{(τ -2 -1) 2 + 4D 1 τ -2 /K}
(X T 0J 0 X 0J 0 /n) | ∞ ≤ |J 0 | 1/2 | (X T 0J 0 X 0J 0 /n) | 2 ≤ (|S * | + |T * |) 1/2
/C min , Theorem 3 now follows from Theorem 1 of Wainwright (2009).

GENERALIZATION OF LEMMA 1

Here, we do not consider the orthogonal and balanced setting anymore and present general conditions ensuring that the irrepresentability conditions (IC) and (IC) 0 are fulfilled by the design matrices X and X 0 involved in the basic approach and our proposal, respectively. For all k ∈ [K], we assume that τ k = τ 0 (n k /n) 1/2 for some τ 0 > 0 and that n

-1/2 k X (k) j 2 ≤ 1 for all (j) ∈ [p].
For X equal to either X or X 0 , this ensures that n -1 X j 2 ≤ max(1, τ -1 ), for each column X j of X .

For any given vector of reference strata ∈ [K] p and any j ∈

[p], set K ,j = {k ∈ [K] : β * k,j = β * j ,j } and K ,j = K ,j \ { j }. Further set, for any k ∈ [K], T ,k = {j ∈ [p] : β * k,j = β * j ,j }, Σ ,k = X (k)T T ,k X (k) T ,k , Π ,k = X (k) T ,k Σ -1 ,k X (k)T T ,k , and
Z (k) ,j = (I n k -Π ,k )X (k) j . Define ω (k) ,j = Σ -1 ,k X (k) T T ,k X (k) j
and Ω

(k) = Σ -1 ,k X (k) T T ,k X (k) S . Introduce the quantities Σ = k∈[K] X (k)T S (I n k -Π ,k )X (k)
S and Ω (k)

,j = Σ -1 X (k)T S Z (k) ,j . Finally set c 1 ( ) = max j∈S c    k∈[K] Ω (k) ,j 1 + k∈[K] τ k l∈[K] Ω (l) Ω (k) ,j 1    c 2 ( ) = max j∈[p] max k∈K ,j    Ω (k) ,j 1 τ k + l =k τ l τ k Ω (l) Ω (k) ,j 1 + ω (k) ,j + Ω (k) Ω (k) ,j 1    c2 ( ) = max j∈[p] max k∈ K ,j    Ω (k) ,j 1 τ k + l =k τ l τ k Ω (l) Ω (k) ,j 1 + ω (k) ,j Ω (k) + Ω (k) ,j 1    . Lemma 4 Let ∈ [K] p be a given vector of reference strata. Assume that Λ min (Σ ,k ) > 0 for k ∈ [K] and Λ min ( Σ ) > 0. Condition (IC ) holds if and only if c 1 ( ) < 1 and c 2 ( ) < 1. Assume that Λ min (Σ * ,k ) > 0 for k ∈ [K] and Λ min ( Σ * ) > 0. Condition (IC) 0 holds if and only if c 1 ( * ) < 1 and c2 ( * ) < 1.
The proof of Lemma 4 follows from the same arguments as those presented in Section 5.1.1 for the proof of Lemma 1, and is omitted. Again, the conditions ensuring that (IC) * and (IC) 0 hold are very similar. This shows that our proposal is able to mimic the optimal version of the basic approach even when the designs are not orthogonal or strata are not balanced.

OTHER RESULTS

Corollary 5 and Corollary 6 consider the two special cases where T * = ∅ and S * = ∅, respectively.

Corollary 5 Assume that the noise variables (ε [K] are independent and identically centered sub-Gaussian variables with parameter σ > 0. Define X = (X (1) T , . . . , X (K) T ) T , the n × p matrix with all the strata pooled together. Set τ k = τ 0 (n k /n) 1/2 > 0 for all k ∈ [K], for some τ 0 > 0. For all j ∈ [p], assume that there exists some β * j ∈ IR such that β * k,j = β * j for all k ∈ [K] and set S * = {j ∈ [p] : β * j = 0}. Further assume that Λ min (X T S * X S * /n) ≥ C min for some C min > 0, and that n

(k) i ) i∈[n k ],k∈
-1/2 k X (k) j 2 ≤ 1 for all (k, j) ∈ [K] × [p].
Finally assume that the three following conditions hold:

( A) k∈[K] τ k > 1, ( C.i.1) c1 := max j / ∈S * (X T S * X S * ) -1 X T S * X j 1 < 1, ( C.ii) c2 := max j∈[p] max k∈[K] τ -1 k (X T S * X S * ) -1 X (k) T S * X (k) j 1 < 1. Now, set γ = (1 -c1 ) ∧ (1 -c2 ) and λ 1 = 2 (1 ∧ τ 0 )γ 2σ 2 log((K + 1)p) n 1/2 , β min = λ 1 |S * | 1/2 C min + 4σC -1/2 min .
Then our proposal run with parameters λ 1 and λ 2,k = τ k λ 1 identifies S * and T * = ∅ with probability at least 1 -4 exp(-c 1 nλ 2 1 ) for some c 1 > 0, as long as min j∈S * |β * j | > β min .

Condition ( C.i.1) is exactly the irrepresentability condition on matrix X, while conditions ( C.ii) and ( A), which are very similar, both simply require that τ 0 is high enough. Moreover, γ = 1 -c1 and 1 ∧ τ 0 = 1 for τ high enough. Therefore, our proposal mimics the lasso run on (X, Y ) provided τ 0 high enough and is optimal, up to log-terms, when the β * k 's are all equal.

Corollary 6 Assume that the noise variables (ε

(k) i ) i∈[n k ],k∈[K]
are independent and identically centered sub-Gaussian variables with parameter σ > 0.

Set τ k = τ 0 (n k /n) 1/2 > 0 for all k ∈ [K], for some τ 0 > 0. For all j ∈ [p], set K * j = {k ∈ [K] : β * k,j = 0} and for all k ∈ [K], set T * k = {j ∈ [p] : β * k,j = 0}. Assume that min k {Λ min (X (k) T T * k X (k) T * k /n k )} ≥ C min for some C min > 0, and that n -1/2 k X (k) j 2 ≤ 1 for all (k, j) ∈ [K] × [p]
. Further assume that the three following conditions hold:

( Ā) ∀j ∈ [p], / ∈K * j τ < 1 + k∈K * j τ k , ( C.i.1) c1 := max j∈[p] max k∈[K] k∈[K] τ k (X (k) T T * k X (k) T * k ) -1 X (k) T T * k X (k) j 1 < 1, ( C.ii) c2 := max k∈[K] max j / ∈T * k (X (k) T T * k X (k) T * k ) -1 X (k) T T * k X (k) j 1 < 1. Now, set γ = (1 -c1 ) ∧ (1 -c2 ), and
λ 1 = 2 (1 ∧ τ 0 )γ 2σ 2 log((K + 1)p) n 1/2 , β min = λ 1 τ 0 |T * | 1/2 C min + 4σC -1/2 min .
Then our proposal run with parameters λ 1 and λ 2,k = τ k λ 1 recovers S * = ∅ and T * = {(k, j) : β * k,j = 0} with probability at least 1-4 exp(-c 1 nλ 2 1 ), for some c 1 > 0, as long as min

(k,j)∈T * |β * k,j | > β min (n/n k ) 1/2 .
Condition ( C.ii) is exactly the union of the irrepresentability conditions for each matrix X (k) , k ∈ [K], while conditions ( C.i.1) and Ā both simply require that τ 0 is small enough. For τ 0 small enough, our proposal then mimics the strategy consisting in performing K lasso on the data (X (k) , y (k) ), k = 1, . . . , K, independently, with a common λ 1 value for each lasso. It is optimal, up to log-terms, when S * = ∅.

Additional empirical results

UNDER THE DESIGNS CONSIDERED IN THE MAIN TEXT

Figure 3 presents additional results regarding the recovery of the set S * 1,P 0 = {j ∈ P 0 : β * 1,j = 0} for the basic approach run with reference strata = (1, . . . , 1) and the recovery of the set S * P 0 = {j ∈ P 0 : β * * j ,j = 0} for the other approaches. Overall, our proposal and the optimal version of the basic approach perform similarly according to this criterion too. In the constant δ * k,j case, all methods perform similarly, and their performance does not depend on either p or d H . In the random δ * k,j case, the performance of each method decreases as p and/or d H increases. This discrepancy with the results obtained in the constant case illustrates that it is harder for the optimal version of the basic approach, our proposal and the clique-based strategy to determine whether the overall effect β * * j ,j of any covariate j is null when the collection of values (β * 1,j , . . . , β * K,j ) varies around zero; keep in mind that in the random case, we have either β * k,j = 1 or β * k,j = 1 ± δ * k,j with δ * k,j ∼ U [K 1/2 /2,2K 1/2 ] so that β * k,j can be negative. Interestingly, it is harder for the basic approach to determine whether β * 1,j is null in this situation too. In addition, the basic approach is generally outperformed by the other approaches in this random δ * k,j case. The clique-based strategy performs well for d H ≥ 6, especially when n k /p is not too small. But it is outperformed by our proposal, and the optimal version of the basic approach, for d H = 3 and p = 500.

UNDER AN ALTERNATIVE SCENARIO

Here, we present additional empirical results obtained under a scenario which should favor the clique-based strategy. We still consider the case where K = 20 and P 0 ⊂ [p], with |P 0 | = 20 but, for each j ∈ [P 0 ], we set β 1,j = • • • = β 5,j = -a, β 6,j = • • • = β 10,j = 0, β 11,j = • • • = β 15,j = a and β 16,j = • • • = β 20,j = 2a, for some a > 0. In other words, for each j ∈ [P 0 ], the effects of the jth covariate across the 20 strata are made of 4 groups of distinct values, which should favor the clique-based strategy. Two values of a were considered, a = K 1/2 and a = K 1/2 /3. Because results were very similar, only those obtained for a = K 1/2 /3 are presented here. Figure 4 presents the predictive performance of each approach. We especially observe that our proposal performs nearly as well as the clique-based strategy in general, and outperforms it for p = 500. It also slightly outperforms the two versions of the basic approach. Figure 5 presents the estimates returned by each approach for one particular simulation in the configuration n k = 100 and p = 100. It can especially be seen that the clique-based strategy does not fuse coefficients sensibly better than the other approaches. Each of the four remaining columns presents the estimates obtained according to one of the four considered approaches: our proposal, the clique-based approach and the basic approach with the reference stratum set to either 1 or 20 for every covariate.

Figure 1 :

 1 Figure 1: Results from the simulation study. The top three panels show the accuracy regarding the recovery of the set T * 1,P 0 = {(k, j) ∈ [K] × P 0 : β * k,j = β * 1,j } for the basic approach and the set T * P 0 = {(k, j) ∈ [K] × P 0 : β * k,j = β *

Figure 2 :

 2 Figure2: Estimation of the K = 8 parameter vectors in the linear regression models describing the association between EGR2 and the p = 44 other transcription factors, at times H0, H1, H6, H12, H24, H48, H72 et H96. Each column corresponds to the estimation obtained according to one of the five considered approaches : the basic approach with the reference stratum set to either H0 or H96 for every covariate, two versions of the generalized fused lasso estimates, one based on cliques and one based on chain graphs, and our proposal.

For

  any s ∈ [|S * |] and t ∈ [|T * k |], denote the sth element of S * by S * s and the tth element of T * k by T * k,t . For all j ∈ S * , further denote by N * j = n|K * j |/K the number of observations in strata contained in K * j = {k ∈ [K] : β * k,j = β * * j ,j }. Now introduce C, the matrix of size |T * | × |S * | made of K blocks C k . Each block is of size |T * k | × |S * |, and the element

Figure 3 :

 3 Figure 3: Accuracy regarding the recovery the set S * 1,P 0 = {j ∈ P 0 : β * 1,j = 0} for the basic approach and of the set S * P 0 = {j ∈ P 0 : β * * j ,j = 0} for the three other approaches. (Left): Constant δ * k,j case. (Right): Random δ * k,j case. Results correspond to averages over 50 replicates in each configuration. Solid line: our proposal. Dotted line: optimal version of the basic approach. Dash-dot line: basic approach. Dashed line: clique-based approach.

Figure 5 :

 5 Figure 5: Estimation of the K = 20 parameter vectors in one particular simulation with n k = 100and p = 100. The first column presents the true values. Each of the four remaining columns presents the estimates obtained according to one of the four considered approaches: our proposal, the clique-based approach and the basic approach with the reference stratum set to either 1 or 20 for every covariate.

  . Sparsity and smoothness via the fused lasso. J. R. Statist. Soc. B, 67, 91-108. VIALLON, V., LAMBERT-LACROIX, S., HOEFLING, H. & PICARD, F. (2016). On the robustness of the generalized fused lasso to prior specifications. Stat Comput, 26, 285-301. VODUC, K. D., CHEANG, M. C., TYLDESLEY, S., GELMON, K., NIELSEN, T.O. & KENNECKE, H. (2010). Breast cancer subtypes and the risk of local and regional relapse. J. Clin. Oncol., 28, 1684-91. WAINWRIGHT, M. (2009). Sharp thresholds for high-dimensional and noisy sparsity recovery using-constrained quadratic programming (lasso). IEEE Trans. Inf. Theory, 55, 2183-202. ZHAO, P. & YU, B. On Model Selection Consistency of Lasso. J. Mach. Learn. Res., 7, 2541-63. ZOU, H. (2006). The Adaptive Lasso and Its Oracle Properties. J. Am. Statist. Assoc., 101, 141829.

  t 2 = S * s for some s ∈ [|S * |], and 0 otherwise. For k ∈ [K], denote by B k,k the diagonal matrix of size |T * k | × |T * k | with tth diagonal term equal to τ 2 k (N * S * s + n k )/(N * s n k ) if T * k,t = S * s for some s ∈ [|S * |] and τ 2 k /n k otherwise. Finally denote by D the diagonal matrix of size |S * | × |S * | with jth diagonal term equal to 1/N * j , and by B the matrix of size |T * | × |T * | made of K 2 blocks, with block (k 1 , k 2 ) equal to B k 1 ,k 2 . By standard algebra, we have

  1/2 . Denote by | M | ∞ the maximum row sum matrix norm of matrix M , and by | M | 2 its spectral norm. Because |

/n, giving overall signalto-noise ratio equal to 1. For each particular choice of n k ∈ {10, 50, 100}, p ∈ {20, 100, 500} and d H ∈ {1

, 3, 6, 9}, and for both the random and constant choice for the δ * k,j 's, we generate 50

Acknowledgement

We are grateful to Adeline Samson and Philippe Rigollet and the anonymous referees and associate editor for their fruitful comments on a preliminary version of this article.

Supplementary material

Supplementary material available at Biometrika online includes two sections. Section 1 presents technical details: the proof of Lemma 1 and a generalized version of Lemma 1, the version of Theorem 1 in the balanced and orthogonal design along with its proof, and two corollaries describing the particular cases where S * = ∅ and T * = ∅. Section 2 presents additional results from our empirical study: accuracies for the recovery of other sets of interest in the settings described heres, and additional results obtained under an alternative settings which should favor the clique-based strategy.

Theorem 3 For all k ∈ [K], assume that the noise variables ε

are independent and identically distributed centered sub-Gaussian variables with parameter σ > 0. Further assuming that (sIC) 0 holds, define

Finally introduce β

min ], and consider the following β-min conditions:

Then, S * and T * are both recovered • with probability superior to 1 -4 exp(-c 1 nλ

1 ), for some c 1 > 0, by the optimal version of the basic approach run with λ 1 = λ

min ;

• with probability superior to 1 -4 exp(-c 1 nλ

(1) 2 1 ), for some c 1 > 0, by our approach run with λ 1 = λ

(1)

min .

Consider the asymptotic setting with K, and possibly p, tending to infinity as n → ∞. Further assume that D * ,0 = D * ,1 = D * . If D * K 1/2 or D * = cK 1/2 for some 0 < c ≤ 1/2, then τ 0 = 1 ensures perfect recovery for signals such that β

), which is optimal up to log-terms. If D * = cK 1/2 for some c > 1/2, we get the same order of magnitude for β

D * , then the regime changes. For K 1/2 D * K the optimal choice is τ 0 = K 1/2 /(2D * ) which only ensures perfect recovery for β

cK for some 0 < c < 1/3, then the result of Theorem S1 becomes almost meaningless: the optimal τ 0 is O(K -1/2 ) which only ensures perfect recovery for β

PROOF OF THEOREM 3

In view of Lemma 1, and because X * J * = X 0J 0 , we simply have to compute Λ min ((X T 0J 0 X 0J 0 /n)) in order to apply Theorem 1 of Wainwright (2009) and establish Theorem 3. To do so, we look for (20, . . . , 20). Dash-dot line: basic approach with = (1, . . . , 1). Dashed line: cliquebased approach.


Jan Dupont

Xavier Orlik
email: xavier.orlik@onera.fr

Romain Ceolato

Thibault Dartigalongue

Spectralon spatial depolarization: towards an intrinsic characterization using a novel phase shift distribution analysis

Keywords:

teaching and research institutions in France or abroad, or from public or private research centers.

Spectralon spatial depolarization: towards an intrinsic characterization using a novel phase shift distribution analysis

Introduction

Identifying the light scattering regime in materials is of interest for various applications, as remote sensing [START_REF] Derode | Robust acoustic time reversal with high-order multiple scattering[END_REF], Lidar measurements [START_REF] Kavaya | Target reflectance measurements for calibration of lidar atmospheric backscatter data[END_REF][START_REF] Kavaya | Polarization effects on hard target calibration of lidar systems[END_REF], biological imaging [START_REF] Abou Nader | Early diagnosis of teeth erosion using polarized laser speckle imaging[END_REF][START_REF] Ceolato | Spectral degree of linear polarization of light from healthy skin and melanoma[END_REF], astrophysics [START_REF] Seager | Photometric light curves and polarization of close-in extrasolar giant planets[END_REF][START_REF] Kervella | The dust disk and companion of the nearby AGB star L 2 Puppis SPHERE/ZIMPOL polarimetric imaging at visible wavelengths[END_REF], and for a better understanding of light matter interactions [START_REF] Sorrentini | Statistical signatures of random media and their correlation to polarization properties[END_REF]. Experimentally, the classical method to measure the linear depolarization produced by a given sample is to illuminate it with a fully linearly polarized light, and to measure the back scattered intensities polarized along the parallel and perpendicular states of polarization. The DOLP is obtained with DOLP = I -I ⊥ I +I ⊥ , and is equal to 0 for a pure depolarizer and close to 1 for a mirror.

Spectralons are reference radiometric samples: their reflectance values are calibrated. Yet, their depolarization rate is rarely fully characterized. This is a major issue as most polarimetric dependency of optoelectronic devices transfer functions must be characterized, especially when using mirrors, diffraction gratings, lasers, etc. Spectralons are used for that purpose, thus there is a need in spectralon polarimetric characterization, in a large range of reflectances (2 -99 %).

In the case of low reflectance Spectralons, the depolarization decrease because of the absorption in the sample: the majority of multiple scattering events are absorbed. Conceptually, the degree of polarization of the back-scattered intensity is expected to reach 1 when the reflectance of the surface is reaching 0%. Measuring such a weak depolarization using the classical method is not trivial as the signal to noise ratio decreases for low reflectances. The polarimetric characterization of high reflectance Spectralon (R = 99%) has been performed through Mueller Matrix measurements [START_REF] Svensen | Mueller matrix measurements and modeling pertaining to Spectralon white reflectance standards[END_REF][START_REF] Sanz | Comprehensive polarimetric analysis of Spectralon white reflectance standard in a wide visible range[END_REF][START_REF] Kildemo | Parametric model of the Mueller matrix of a Spectralon white reflectance standard deduced by polar decomposition techniques[END_REF], resulting in a quasi perfect depolarizing Mueller Matrix, for small illumination and detection angles. However, an extrinsic spatial depolarization is often generated to overcome the speckle noise issue: the electromagnetic field is (temporally or spatially) averaged during measurements. That decrease in degree of polarization due to an incoherent spatial summation of various contibutions to the electromagnetic field has been for example observed in [START_REF] Li | Degree of polarization in laser speckles from turbid media: Implications in tissue optics[END_REF][START_REF] Ghabbach | Depolarization and enpolarization DOP histograms measured for surface and bulk speckle patterns[END_REF], and studied in the case of surface [START_REF] Zerrad | Gradual loss of polarization in light scattered from rough surfaces: Electromagnetic prediction[END_REF] and bulk [START_REF] Zerrad | Spatial depolarization of light from the bulks: electromagnetic prediction[END_REF] scattering. The formalism describing the Mueller matrix dependency with the spatial coherence of the fields has been studied in [START_REF] Ossikovski | General formalism for partial spatial coherence in reflection Mueller matrix polarimetry[END_REF].

In this paper, we propose to perform polarimetric measurements in the spatially resolved speckle fields, i.e. close to the limiting case of a total spatial coherence. Thus, the spatial depolarization is limited, and we approach intrinsic polarimetric parameters characterizing the samples, that does not depend on imaging processes or other external parameters. To perform these measurements, we use the SOPAFP method [START_REF] Dupont | Polarization analysis of speckle field below its transverse correlation width : application to surface and bulk scattering[END_REF][START_REF] Dupont | Polarized vortices in optical speckle field: observation of rare polarization singularities[END_REF] which increases the polarimetric accuracy, especially at places of destructive interferences, and in case of low reflectance samples. As a result, we are able to fully characterize the polarization state for every pixel of the speckle image. Every pixel exhibits a linear combination of a fully depolarized light, and a fully polarized light. The first term includes experimental noises and local depolarization due to eventual fields integration on the pixel surface, if the speckle is not enough spatially resolved. In order to study the depolarization produced by the samples, we spatially average the polarized fields on increasing areas, which globally lowers the DOP. This phenomenon is called spatial depolarization and leads to a first data set of DOP as a function of spatial averaging, the asymptotic values being the apparent DOP values that one would measure by integrating the whole fields on a photodiode surface, i.e. close to the limiting case of total spatial incoherence.

In order to increase the accuracy of the spatial depolarization determination, we come back to the measurements obtained in the case of a quasi total coherence (i.e. in spatially resolved speckle fields), then the statistics of the relative phase shifts between two orthogonal components of the electromagnetic field are studied. We perform the projection of all the Stokes vectors on the plane waves polarized at ±45 • relatively to the mean scattered state. As the Lambertian samples does not exhibit spatial variations of diattenuation, the mean intensities of the plane waves are equal, and the SOP dispersion around the mean state is only driven by the relative phase shift probability density function (PDF). The relative phase shifts PDF appears to be close to circular distributions [START_REF] Gatto | The generalized Von Mises distribution[END_REF][START_REF] Mccullagh | Möbius transformation and Cauchy parameter estimation[END_REF][START_REF] Elderton | Systems of Frequency Curves[END_REF][START_REF] Hall | The approximation of symmetric X-ray peaks by Pearson type VII distributions[END_REF]. The bandwidth of the distribution, depending on a concentration parameter ρ, is demonstrated to be equal to the spatially integrated DOP. For Lambertian samples, the concentration of the phase distribution ρ appears to be only depending on the reflectance R of the samples. Thus, an analytic equation is given, allowing to estimate the SOP dispersion, or the monochromatic spatially integrated DOP (if the fields are incoherently spatially averaged), as a function of the reflectance of the Lambertian samples.

Method and formalism

Any State Of Polarization (SOP) can be described by a sum of two plane waves, orthogonal to each other, with amplitudes E x and E y , and φ their relative phase shift. If we consider a polarized wave propagating along the z axis, its SOP is described by the following Stokes vector [START_REF] Goldstein | Polarized Light[END_REF]:

S =            I = E 2 x + E 2 y + C 2 Q = 2 E x E y cos (φ) U = E 2 x -E 2 y V = 2 E x E y sin (φ) (1)
With I the total intensity (polarized or not), Q the intensity polarized along the horizontal or vertical, U the intensity polarized in the ±45 • direction and V the intensity polarized along a right or left handed circular SOP; E x the amplitude of the first plane wave linearly polarized along the x axis, E y the amplitude of the second plane wave linearly polarized along the y axis (with x and y corresponding respectively to ±45 • direction), φ their relative phase shift, and C 2 the constant intensity corresponding to the depolarized fields and experimental noises [START_REF] Dupont | Speckle fields polarimetry: statistical analysis and polarization singularities measurements[END_REF]. The corresponding local Degree Of Polarization (DOP) is:

DOP = Q 2 + U 2 + V 2 I = 1 - C 2 I (2)
In order to perform the polarimetry of the scattered fields, the State Of Polarization Analysis by Full Projection on the Poincaré sphere (SOPAFP) method [START_REF] Dupont | Polarization analysis of speckle field below its transverse correlation width : application to surface and bulk scattering[END_REF][START_REF] Dupont | Polarized vortices in optical speckle field: observation of rare polarization singularities[END_REF] is used, which consists in the non linear regression of the intensity recorded during n polarimetric projections. These multiple projections (300 for the measurements presented in this paper) increase the polarimetric accuracy, especially at places of destructive interference, where the signal to noise ratio drastically decreases.

We define the two intensities I x = |E x | 2 and I y = |E y | 2 , the intensities of the two plane waves respectively polarized along the x and y directions, and the following intensity ratios:

R x = < I x > < I > R y = < I y > < I > R c = < C 2 > < I > (3)
With <> indicating the spatial averaging. R x describes the mean fraction of the total intensity polarized along the x direction; R y the mean fraction of the total intensity polarized along the y direction; and R c the mean fraction of totally depolarized intensity.

The spatially integrated degree of polarization DOP, not equivalent to the mean local degree of polarization < DOP >, is the degree of polarization that one would measure after a spatial integration of all the fields incoming in the CCD, e.g. by replacing it with a photo-diode. The spatially integrated DOP can be computed with:

DOP = < Q > 2 + < U > 2 + < V > 2 < I > (4)
With I, Q, U and V described in Eq. ( 1). In order to measure the spatial correlation between the two intensities I x and I y , we use the sample Pearson correlation coefficient (r), defined by:

r xy = m n (x mn -x)(y mn -ȳ) ( m n (x mn -x) 2 )( m n (y mn -ȳ) 2 ) (5)
With m, n the pixels coordinates, and -1 ≤ r xy ≤ 1. If r xy = -1 or r xy = 1 the two intensities exhibit a negative or positive spatial correlation, if r xy = 0 the two intensities are spatially uncorrelated.

Finally, we describe the relative phase shift PDF p(φ) with different circular distributions. We made a first attempt with a Von Mises distribution V M(φ), approximation of a circular normal distribution [START_REF] Gatto | The generalized Von Mises distribution[END_REF]:

V M(φ) = 1 2π exp(κ.cos(φ -φ 0 )) I 0 (κ) , 0 ≤ φ < 2π (6)
With κ a parameter describing the concentration of the Von Mises distribution, I n (x) the n th order modified Bessel function and φ 0 a location parameter. We also tried to describe the relative phase shift distribution p(φ) with a circular Cauchy (or Lorentzian) CC(φ) distribution [START_REF] Mccullagh | Möbius transformation and Cauchy parameter estimation[END_REF]:

CC(φ) = 1 2π 1 -ρ 2 1 + ρ 2 -2ρcos(φ -φ 0 ) , 0 ≤ φ < 2π (7)
With ρ the concentration of the circular Cauchy distribution. In order to get a better accuracy, we finally tried the following circular convolution:

CV = V M ⊗ CC (8)
Where CV denotes circular Voigt, analogy of a Voigt profile, well defined as the convolution of a classical Gaussian and Lorentzian (also called Cauchy) profile [START_REF] Elderton | Systems of Frequency Curves[END_REF][START_REF] Hall | The approximation of symmetric X-ray peaks by Pearson type VII distributions[END_REF].

Experimental setup

We use the experimental set up presented in Fig. 1. The wavelength is 532 nm. The illumination beam is collimated, with a diameter of roughly 1 cm, and the incident SOP is set to the vertical by a Glan-Taylor polarizer (P lin ). The illumination beam is scattered by the sample, the polarimetric projections of the scattered field are performed by two Nematic Liquid Crystals (N LC) and a linear grid polarizer (P lin ). The resulting intensity is imaged by a lens, which entrance pupil is reduced by a pinhole (P), on a 1000×1000 pixels CCD. The distance between the sample and the entrance pupil is roughly 20 cm, the measurement is performed in the far field region. The Fig. 1. Schema of the experimental setup. The source is a SLM laser @532nm. The illumination SOP is set to the vertical by P lin . The samples are Lambertian Spectralons, the polarimetric projections of the scattered fields are performed by two Nematic Liquid Crystals (NLC) and a vertical polarizer. The resulting intensity is imaged by a lens, in which the imaging pupil is reduced using a diaphragm P, and a CCD camera.

illumination is perpendicular to the sample surface, and the detection is centered roughly at 45 • respectively to the illumination beam. We perform measurements on Lambertian samples with calibrated reflectances. Thus, the surface luminance is approximately invariant in function of the detection angle. The designation of the samples are SRS -X X -010, with X X the theoretical reflectance value, from 02% to 99%. We also use the colored sample SCS -RD -010. The reflectances are measured using a PerkinElmer Spectrophotometer @ 532 nm, the results are resumed in Table 1 for each sample, here we describe the samples through their reflectance values.

Table 1. Diffuse reflectance R of each sample measured using a PerkinElmer Spectrophotometer @ 532 nm.

Sample SRS-02 SCS-RD SRS-20 SRS-50 SRS-75 SRS-99 R (%)

1.2 6.9 18.4 49.9 77.4 98.7

In order to reduce the sources of depolarization (i.e. to reduce the C 2 value), some precautions are taken: the laser source is Single Longitudinal Mode (SLM), in order to avoid spectral depolarization, and we reduce the diameter of the imaging pupil, in order to reduce the pixel spatial depolarization. As the diameter of the imaging pupil is reduced, the size of the imaging Point Spread Function is increased by the diffraction. Thus, the number of pixels per speckle grain in the image plane is increased and we have a correct spatial sampling of the diffraction pattern. Indeed, the quality of the spatial sampling impacts the amount of spatial depolarization [START_REF] Zerrad | Gradual loss of polarization in light scattered from rough surfaces: Electromagnetic prediction[END_REF][START_REF] Zerrad | Spatial depolarization of light from the bulks: electromagnetic prediction[END_REF].

Thus, the contributions remaining in the C 2 value are the intensity variations uncorrelated with polarimetric projections, mostly experimental noises as photon noise, non linearity of the pixels quantum efficiency, dark current, etc. As the experimental noises can lead to higher or lower DOP values, the C 2 value is not restricted to be positive, i.e. C is a complex number.

Results

SOP statistics: spatial depolarization measurements

We perform polarimetric measurement using the SOPAFP method [START_REF] Dupont | Polarization analysis of speckle field below its transverse correlation width : application to surface and bulk scattering[END_REF][START_REF] Dupont | Polarized vortices in optical speckle field: observation of rare polarization singularities[END_REF] on reference samples, with calibrated reflectances. As the scattered fields are temporally stable, we neglect the impact of the time of exposure on polarimetry, and we adjust the integration time of the CCD in order to avoid the saturation, the illumination beam being unchanged for all samples. The intensities I and C 2 , the relative phase shifts φ and the SOP probability density in the Poincaré sphere measured in the fields scattered by the samples with R = 1%, R = 7% and R = 99% are displayed in Fig. 2, for a 50×50 pixels area. By looking at the relative phase shifts φ, Fig. 2(c) shows that the sample R = 1% scatters a field in which phase shifts are centered around π, while no predominant value can be distinguished in the field scattered by the sample R = 99% in Fig. 2(k). Moreover, the SOP are spreading on the Poincaré sphere surface when the sample reflectance increases: a majority of the scattered SOP are centered around the illumination SOP for the sample R = 1%, Fig. 2(d), while there is no predominant SOP scattered by the sample R = 99%, Fig. 2(l). With the plane waves decomposition described by Eq. ( 1), and because the illumination SOP is set to the vertical (i.e.

I (n.c.) C 2 (n.c.) φ (rad) PDF(S) (arb.) R = 1% (a) (b) (c) (d) R = 7% (e) (f) (g) (h) R = 99% (i) (j) (k) (l)
x and y are respectively at ±45 • relatively to the mean scattered state), for each sample we obtain exponentially decreasing Probability Density Functions (PDF) for p(I x ) and p(I y ), with equal mean intensity levels: R x = R y . Thus for all the samples, I x and I y are two independent speckle fields [START_REF] Goodman | Speckle Phenomena in Optics[END_REF], with equal mean intensities. Reciprocally, in order to modelize the polarized subjective speckle produced by any kind of Lambertian scatterer (from a simple to a multiple scatterer), a simple way is to sum two independent speckles with the same mean intensity, orthogonally polarized and with partially correlated phase planes, the scattering regime depending on that partial correlation [START_REF] Dupont | Simulation of polarized optical speckle fields: effects of the observation scale on polarimetry[END_REF].

The intensity spatial correlation r xy between I x and I y is equal to r xy = 0.72 for the sample R = 1%, and r xy = 0.14 for the sample R = 7%. For samples with higher reflectances, the spatial correlation falls roughly around 0 and is spatially variant. Thus this parameter is not sufficiently sensitive to fully describe the polarimetric behavior of the samples on the whole reflectance range 1% -99%. In classical polarimetry, usually the lens aperture is large relatively to the wavelength, and the effect of diffraction, the speckle pattern, is spatially summed on the pixels surface. This summation leads to a measurement of the DOP of several contributions to the field, spatially integrated. For our measurements, the surface of a single speckle grain is roughly 15 × 15 pixels. The corresponding f-number N can be estimated through the speckle transverse width l x in the approximation of a diffraction limited lens: N l x /(2.44λ). For 4.4 µm square pixels, and a speckle transverse width of 15 pixels, the f-number is approximated by N 51.

For all samples, we measure 3% < R c < 7%, thus the mean local DOP value measured in the speckle fields remain relatively high. We spatially sum the measured SOP on larger areas, from 2 × 2 to 300 × 300 pixels windows, leading to a deacreasing of the apparent spatially integrated Degree Of Polarization (DOP). The DOP variations in function of the size of spatial averaging is displayed in Fig. 3. The DOP isn't monotonically decreasing in function of the size of the spatial summation, it depends on the spatial gradient of SOP variations and on the C 2 intensity, which are fluctuating at the speckle grain scale. The fact that the DOP is globally decreasing with the spatial integration is called "spatial depolarization" [START_REF] Li | Degree of polarization in laser speckles from turbid media: Implications in tissue optics[END_REF][START_REF] Zerrad | Gradual loss of polarization in light scattered from rough surfaces: Electromagnetic prediction[END_REF][START_REF] Zerrad | Spatial depolarization of light from the bulks: electromagnetic prediction[END_REF]. The asymptotic values (for a spatial averaging of 300 × 300 pixels) are indicated in Fig. 3, and are corresponding to the values measured using a large pupil aperture relatively to the wavelength. For example, in the case of the same sample R = 99%, we obtain a good agreement with the classical measurement method: the spatially integrated DOP (here equal to 0.04) is measured "lower than 0.03" in [START_REF] Ceolato | Spectral degree of linear polarization of light from healthy skin and melanoma[END_REF].

By looking at Fig. 3 (obtained with 15 × 15 pixels per speckle grains) one can notice that the DOP value is stable with at least a 100 × 100 pixels window size, which would corresponds roughly to the integration of 50 speckle grains on each pixel.

The spatially integrated DOPs values DOP are here obtained for one incident SOP (vertical), and for fixed illumination and detection angles relatively to the sample normal (resp. 0 • , 45 • ). The DOP of R99 has been widely studied in function of these angles and its dependency with the illumination SOP through Mueller polarimetry, e.g. in [START_REF] Kildemo | Parametric model of the Mueller matrix of a Spectralon white reflectance standard deduced by polar decomposition techniques[END_REF]. Here, we propose to return to the measurements obtained in the spatially resolved speckle fields (i.e. in the case of quasi total coherence, leading to the local SOP and DOP values), and to characterize, for the experimental setup used here, the process leading to that apparent DOP decreasing, and its dependency with the sample reflectance.

Relative phase shifts statistics: characterization of the intrinsic spatial depolarization

Figure 4 displays the PDFs of the relative phase shifts weighted by their corresponding intensity values p(φ) for all samples, measured in their respective speckle fields, and the one measured in a field scattered by a mirror. As a result, distributions centered around the relative phase shift π are observed, which correspond, when I x = I y , to the vertical incident polarization state: )). Moreover, one can see that when the reflectance increase, p(φ) becomes broader. The evolution of the measured standard deviation σ of the relative phase shifts in function of the sample reflectance is displayed in Fig. 5 (black crosses). This evolution can be fitted by the following equation:

S in = [1 -1 0 0] (see Eq. ( 1
σ(R) = σ M I N + τ σ τ σ + 1/R * (9)
With σ the standard deviation of the relative phase shifts φ, R * = R 100-R , R * ∈ [0 ; +∞[, σ M I N an offset parameter (when the reflectance reaches 0%), and τ σ a scale parameter. By regression of Eq. ( 9) on the experimental results, we obtain σ M I N = 0.77 rad, τ σ = 21.48. The adjusted curve is plotted in Fig. 5 (red). The theoretical maximum value for the standard deviation σ is the one of a uniform (U) PDF in the interval [0 ; 2π]:

σ(U) = 2π √ 12
1.81 rad (black dashed line). Fig. 4. Representation of the measured relative phase shifts probability density functions, weighted by the intensity value, for the samples with various reflectances and a mirror as a reference. The PDFs integral are normalized to 1. One can see that when the reflectance increases, the phase distribution tends to be uniform. Thus, the standard deviation of the relative phase shifts σ is depending on the sample reflectance, i.e. the more the light has undergone multiple scattering events, the more σ increases following Eq. ( 9), and almost reaches its theoretical value: max(σ) = 1.77 σ(U). One can notice that the offset value (σ M I N = 0.77 rad) is representative of a constant phase shift dispersion, whatever the sample reflectance value. The offset value σ M I N is roughly 40 times higher than the phase dispersion in the field scattered by the mirror (σ mirror = 0.02 rad). However, it is known that the polarimetric behavior of a sample which reflectance is reaching 0% should be similar to the one of a quasi perfect surface scatterer, quasi all multiple scattering events being absorbed due to its low reflectance value. Thus, the standard deviation of the relative phase shifts appears to be a biased interval estimator for the phase dispersion quantification.

These distributions appears to be more accurately described by circular distributions, we adjust various distributions on the experimental measurements. The mean error between adjusted distribution y(φ) and the experimental data p(φ) is computed with:

err = φ |y(φ) -p(φ)| φ p(φ) (10)
We adjust circular distributions by regression of Eqs. ( 6), ( 7) and ( 8) on the p(φ) distributions, for each samples. In the case of the sample R = 1%, the error for the p(φ) description using a Von Mises distribution V M(φ) (circular normal distribution) is 13%, the correlation coefficient is R 2 = 0.994. Using a circular Cauchy CC(φ) distribution, the mean error is 23%, the correlation coefficient is R 2 = 0.980. Finally, we tried the convolution of the Von Mises with the circular Cauchy distribution, that we called CV(φ), Circular Voigt, the mean error is 6%, the correlation coefficient is R 2 = 0.999. For higher reflectance samples, we measure quasi equal ρ parameters and errors using a CC and a CV distribution, denoting that these two distributions are quasi equivalent and that the circular Cauchy distribution is sufficient for the p(φ) description for samples with R ≥ 7%. Thus, for the sample R = 1% for which the SNR decreases, we present the ρ parameter measured using a CV distribution, and for higher reflectance samples, we use a CC distribution. The fluctuations of the ρ parameter in function of the reflectance are plotted in Fig. 6 (black crosses), with their mean errors for the p(φ) description (blue dots). Finally, we plot the best adjusted fluctuation (red dashed line) computed with Eq. [START_REF] Kildemo | Parametric model of the Mueller matrix of a Spectralon white reflectance standard deduced by polar decomposition techniques[END_REF]:

ρ(R) = ρ M I N + τ ρ τ ρ + R * (11)
With R * = R 100-R , R * ∈ [0 ; +∞[. One could use Eq. ( 11) in order to estimate the relative phase dispersion produced by any sample corresponding to that type of scatterer, for a vertical incident SOP and normal incidence / 45 • detection, only given its reflectance value. By regression, we measure ρ M I N = 0.05 and τ ρ = 0.05. When the reflectance is reaching 0%, the relative phase concentration is reaching 1.05 (i.e. close to the Dirac function). Moreover, when the reflectance is equal to 100%, the relative phase concentration is reaching 0.05 (close to a uniform distribution). These results have to be compared with the measurements of the phase dispersion using the standard deviation (Fig. 5): using that interval estimator, σ M I N = 0.77 rad for R = 0% (see Fig. 5). Thus, ρ appears to describe the intrinsic depolarization of the sample, the extrinsic part (experimental noises and artifacts) being extracted by the convolution with a VM distribution, especially in the case of low reflectance samples.

Indeed, we can draw analogy with molecular spectroscopy, where Gaussian profile generally describes inhomogeneous frequency broadening, while the Lorentzian profile describes homogeneous broadening, their convolution leading to the well defined Voigt profile [START_REF] Elderton | Systems of Frequency Curves[END_REF][START_REF] Hall | The approximation of symmetric X-ray peaks by Pearson type VII distributions[END_REF]. Here, the "homogeneous broadening" of relative phase shifts (CC distribution) would only be depending on the sample. However, the "inhomogeneous broadening" of relative phase shifts (VM distribution) would be depending on experimental noises and artifacts (angular diversity due to the incident beam divergence, angular diversity in the detection line, decreasing of the SNR for low reflectances, ...). Thus, the inhomogeneous broadening (VM distribution) is expected to be only dependent on the experimental set-up. Moreover, in [START_REF] Ossikovski | General formalism for partial spatial coherence in reflection Mueller matrix polarimetry[END_REF] the polarimetric formalism dependant on the illumination / detection spatial coherence is found to be the convolution of the material response and an instrumental complex phase function.

Thus by studying only the homogeneous part (CC distribution) of the circular Voigt profile (convolution of CC with VM), we believe to study only the intrinsic part of the relative phase shifts dispersion, i.e. only due to the scattering by the sample.

Besides, the DOP value can be computed with Eq. ( 2). Back to Eq. ( 1) with C = 0, E x = E y = 1/ √ 2 (as the diattenuation of the samples can be neglected), it comes I(φ) = 1, Q(φ) = cos(φ), U(φ) = 0 and V(φ) = sin(φ). By integration on the entire field (i.e. all the φ values), one obtain:

Q = ∫ 2π 0 cos (φ) p(φ) dφ V = ∫ 2π 0 sin (φ) p(φ) dφ (12)
By substituting p(φ) by Eq. ( 7), i.e. p(φ) = CC(φ) in Eq. ( 12), one obtain V = 0 and Q = ρ. Thus, Eq. ( 2) can be rewritten as: Which indicates that the concentration of the relative phase weighted by its corresponding intensity, described by a circular Cauchy distribution is directly equal for that type of material to the DOP value. We plot in Fig. 7 the DOP (the asymptotic DOP values in Fig. 3) of each sample (black crosses), and the measured relative phase concentration ρ (red crosses). One can observe a good agreement between these two parameters. Equation (11) becomes:

DOP = ρ (13)
ρ(R) = DOP(R) = DOP M I N + τ dop τ dop + R * (14)
With DOP M I N = 0.05 = ρ M I N , τ dop = 0.05 = τ ρ . Thus, for that type of material and with our experimental setup, we have shown that the DOP is solely dependent on the sample reflectance. The offset value in Eq. ( 14) is relatively small: DOP M I N = 0.05, which is roughly equivalent to the instrument uncertainties: R c = 7% for the sample R = 99%. Moreover, for a theoretical sample with a reflectance reaching 0%, DOP M AX is reaching 1.05, which is close to 1 given the instrument uncertainties: the mean local DOP value measured in the field reflected by a mirror is < DOP >= 1.02, with a standard deviation σ(DOP) = 0.03. Similarly, for a vertical illumination SOP and similar illumination and detection angles, one could use Eq. ( 14) in order to estimate the spatially integrated DOP of the fields scattered by any equivalent Spectralon, only given its reflectance value.

Additional studies could be done to analyze the potential effects producing the homogeneous and inhomogeneous contributions to the relative phase shifts dispersion, in case of low reflectance samples: e.g. impact of the angle diversity on the sample interface and detection line, effect of the size / shape diversity of the scatterers relatively to the wavelength, size of the pixel relatively to the speckle grain size or the absorption and scattering coefficients µ a and µ s .

Conclusion

In this paper, experimental polarimetric measurements on reference radiometric samples have been performed. The method used allows to characterize the scattered fields with a high spatial resolution, which limits the spatial depolarization effect at the pixel scale.

The relative phase shifts dispersion produced by the scattering events in that type of samples can be described by circular Cauchy distributions. We proposed the characterization of the samples through their relative phase shifts statistics, with the concentration parameter ρ of circular Cauchy distributions, which appears to be a sensible intrinsic parameter. We saw that this ρ value is decreasing with the sample reflectance, thus ρ is decreasing with the number of multiple scattering events in the material. By regression of the ρ fluctuations in function of the reflectance, the amplitude of phase dispersion produced by multiple scattering in that type of material has been characterized. The ρ measurement can be useful for material distinction, or for defects or inhomogeneities characterization, via the detection of scattering regime variations.

Finally, this ρ parameter appears to be directly equal to the spatially integrated DOP. Thus, the equation allowing to compute the monochromatic DOP value of the entire field scattered by a same type of material with a known reflectance has been calibrated.

Fig. 2 .

 2 Fig. 2. First row: R = 1%, second row: R = 7%, third row: R = 99%. (a, e, i) Intensity I scattered by the samples, expressed in numerical counts. (b, f, j) Totally depolarized intensity and experimental noises C 2 , expressed in numerical counts. (c, g, k) Relative phase shifts φ between the orthogonal plane waves polarized along the x and y directions. (d, h, l) Probability density functions of the measured SOPs represented in the Poincaré sphere.

Fig. 3 .

 3 Fig. 3. Evolution of the measured spatially integrated DOP in function of the size of the window of spatial summation, for the 6 different samples with various reflectances. The asymptotic values only depend on the sample reflectance, and correspond to the values that one would measure with a large pupil aperture relatively to the wavelength.

Fig. 5 .

 5 Fig. 5. Evolution of the measured standard deviation σ of the relative phase shifts φ in function of the reflectance R of the samples (black crosses). Regression of Eq. (9) on experimental measurements (red). Standard deviation of an uniform distribution on [0 ; 2π] (black dashed line).

Fig. 6 .

 6 Fig. 6. Measured ρ parameters (cross, black) of respectively the best circular Voigt (R < 7%) and circular Cauchy (R ≥ 7%) distributions describing the p(φ) distributions, for each sample. (red) Regression of the ρ fluctuations with Eq. (11). Dots: mean error between the p(φ) distribution and respectively the best circular Voigt (R < 7%) and circular Cauchy (R ≥ 7%) distributions.

Fig. 7 .

 7 Fig.7. Measured ρ values for each sample, using a circular Voigt profile for R < 7% and a circular Cauchy profile for R ≥ 7% (red crosses). For the sample R = 1%, we display the ρ measurement using also a circular Cauchy distribution (blue cross). Measured spatially integrated DOP (black crosses). DOP and ρ best adjusted Eq. (11) (dashed lines, resp. black and red). Mean level of local depolarization and experimental noises R c (green dots). One can see a good agreement between ρ and DOP.

Acknowledgments

We would like to thank the region Midi-Pyrénées for providing us financial support for this work.

Funding

Office National d'Etudes et Recherches Aérospatiales (ONERA); Région Midi-Pyrénées, Project APRTCN.