C Soize
email: christian.soize@univ-paris-est.fr

Optimal partition in terms of independent random vectors of any non-Gaussian vector defined by a set of realizations

Keywords: partition in independent random vectors, high dimension, non-Gaussian AMS subject classifications. 62H10, 62H12, 62H20, 60E10, 60E05, 60G35, 60G60 Euclidean vector space of dimension N

We propose a fast algorithm for constructing an optimal partition, in terms of mutually independent random vectors, of the components of a non-Gaussian random vector that is only defined by a given set of realizations. The method proposed and its objective are different from the independent component analysis (ICA) that was introduced to extract independent source signals from a linear mixture of independent stochastic processes. The algorithm that is proposed is based on the use of the mutual entropy from information theory and on the use of graph theory for constructing an optimal partition. The method has especially been developed for random vectors in high dimension and for which the number of realizations that constitute the data set can be small. The proposed algorithm allows for improving the identification of any stochastic model of a non-Gaussian random vector in high dimension for which a data set is given. Instead of directly constructing a unique stochastic model for which its stochastic dimension, which is identified by solving a statistical inverse problem, can be large, the proposed preprocessing of the data set allows for constructing several mutually independent stochastic models with smaller stochastic dimensions. Consequently, such a method allows for decreasing the cost of the identification and/or to make possible an identification for a case that is a priori in high dimension and that could not be identified through a direct and global approach. The algorithm is completely defined in the paper and can easily be implemented.

M +

n : set of all the positive-definite (n × n) real matrices. (Θ, T , P): probability space. L 2 (Θ, R n ): space of all the R n -valued second-order random variables on (Θ, T , P).

x T : transpose of the column matrix representing a vector x in R N .

[x] kj : entry of matrix x : Euclidean norm of vector x. δ ij : Kronecker's symbol such that δ ij = 0 if i = j and = 1 if i = j. pdf: probability density function.

Introduction.

In this introduction, we explain the objective of the paper, we give some explanations concerning the utility of the statistical tool that is proposed, we present a very brief state of the art, and finally, we give the organization of the paper and we describe the methodology proposed.

(i) Objective. The objective of this paper is to propose an algorithm for constructing an optimal partition of the components of a non-Gaussian second-order random vector H ν = (H ν 1 , . . . H ν n ) in terms of m mutually independent random vectors Y ν,1 , . . . , Y ν,m with values in R µ 1 , . . . , R µm , such that H ν = (Y ν,1 , . . . , Y ν,m ), where m must be identified as the largest possible integer and where µ 1 , . . . , µ m must be identified as the smallest integers greater than or equal to one, such that µ 1 + . . . + µ m = n. For all j = 1, . . . , m, the components Y ν,j 1 , . . . , Y ν,j µ j of the R µ j -random vector Y ν,j are mutually dependent. For such a construction, it is assumed that the solely available information is made up of a given data set made up of ν "experimental" realizations η exp,1 , . . . , η exp,ν of the random vector H ν . In the method proposed, a probability density function η → p H ν (η) on R n will be constructed by the Gaussian kernel estimation method from the set of the ν realizations. Then, the random vector H ν will be defined as the random vector for which its probability distribution will be p H ν (η) dη. This is the reason why the random vector is denoted as H ν (and not H) for indicating its dependence on ν). We are interested in the challenging case for which n can be big (high dimension) and for which the number ν of realizations can be small (the statistical estimators can then not be well converged). Note that the "experimental" realizations can come from experimental measurements or come from numerical simulations. Such an objective requires a priori (i) to construct an adapted criterion for testing the mutual independence of random vectors Y ν,1 , . . . , Y ν,m in a given approximation (the exact independence will not be reached due to the fact that the number ν of "experimental" realizations is finite, and can be small enough), and (ii) to develop an algorithm that is faster than the one that would be based on the exploration all the possible partitions of the n components of H ν (without repetition of the indices). The number of possible partitions that should be constructed and tested with the adapted criterion would then be n-1 j=1 C j n with C j n = ((nj + 1) × . . . × n)/j! , which is really too big for a high value of n. In this paper, we propose a statistical tool with a fast algorithm for solving this problem in high dimension and for which the value of ν can be small enough.

(ii) What would be the utility of the statistical tool proposed? For many years, probabilistic and statistical methodologies and tools have been developed for uncertainty quantification in computational sciences and engineering. In particular advanced methodologies and algorithms are necessary for modeling, solving, and analyzing problems in high dimension. A recent state of the art concerning all these methodologies can be found in the Handbook of Uncertainty Quantification [START_REF]Handbook for Uncertainty Quantification[END_REF]. The statistical tool presented in this paper, which consists in constructing an optimal partition of a non-Gaussian random vector in several mutually independent random vectors, each one having a smaller dimension than the original random vector, can be used for decreasing the numerical cost of many existing methodologies devoted to the construction of stochastic models and their identification by solving statistical inverse problems, both for the uncertain parameters of the computational models and the stochastic representation of their random responses. Let us consider, for instance, the case of the polynomial chaos expansion (PCE) of random vectors, stochastic processes, and random fields, which has been pioneered by Ghanem and Spanos [START_REF] Ghanem | Stochastic Finite Elements: A spectral Approach[END_REF] in the community of computational sciences and engineering, on the basis of the mathematical works performed by Wiener [START_REF] Wiener | The homogeneous chaos[END_REF] and Cameron and Martin [START_REF] Cameron | The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals[END_REF]. This initial work has then been extended by Ghanem and his co-authors and by many researchers in the community of uncertainty quantification, concerning the development of many methodologies and tools, and for which numerous applications have been performed in many fields of science and engineering with great success. In the area of the methodologies, we refer, for instance, to [START_REF] Arnst | Identification of Bayesian posteriors for coefficients of chaos expansion[END_REF][START_REF] Das | Asymptotic sampling distribution for polynomial chaos representation of data : A maximum-entropy and fisher information approach[END_REF][START_REF] Debusschere | Numerical challenges in the use of polynomial chaos representations for stochastic processes[END_REF][START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF][START_REF] Doostan | Stochastic model reduction for chaos representations[END_REF][START_REF] Ghanem | Numerical solution of spectral stochastic finite element systems[END_REF][START_REF] Ghanem | A probability construction of model validation[END_REF][START_REF] Le | Spectral Methods for Uncertainty Quantification with Applications to Computational Fluid Dynamics[END_REF][START_REF] Najm | Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics[END_REF][START_REF] Nouy | Random fields representations for stochastic elliptic boundary value problems and statistical inverse problems[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF] for the spectral approaches of linear and nonlinear stochastic boundary value problems and some associated statistical inverse problems; to [START_REF] Ernst | On the convergence of generalized polynomial chaos expansions[END_REF][START_REF] Lucor | Generalized polynomial chaos and random oscillators[END_REF][START_REF] Soize | Physical systems with random uncertainties : Chaos representation with arbitrary probability measure[END_REF][START_REF] Wan | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF][START_REF] Xiu | Wiener-Askey polynomial chaos for stochastic differential equations[END_REF] for the generalized chaos expansions; to [START_REF] Arnst | Hybrid sampling/spectral method for solving stochastic coupled problems[END_REF][START_REF] Arnst | Reduced chaos expansions with random coefficients in reduced-dimensional stochastic modeling of coupled problems[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Ghanem | Stochastic convergence acceleration through basis enrichment of polynomial chaos expansions[END_REF][START_REF] Keshavarzzadeh | Convergence acceleration of polynomial chaos solutions via sequence transformation[END_REF][START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems[END_REF] for works devoted to the dimension reduction and to the acceleration of stochastic convergence of the PCE. If we were interested in the identification of the PCE for the random vector H ν with values in R n , the use of the statistical tool that is proposed in this paper would allow for obtaining a significant gain if the optimal partition was such that 1 ≪ m ≤ n. We briefly explain hereinafter the reason why. For N d ≥ 1 (N d will be the maximum degree of the polynomials) and for 1 ≤ N g ≤ n (N g will be the length of the random germ for the polynomial chaos), let α = (α 1 , . . . , α Ng ) ∈ N Ng be the multi-index such that 0 ≤ α 1 + . . . + α Ng ≤ N d . The first 1 + K multi-indices are denoted by α (0) , . . . , α (K) in which α (0) is the multi-index (0, . . . , 0), and where K = (N g + N d )!/N g ! N d ! -1. For N g fixed, let {Ψ α (ξ) , α ∈ (α 1 , . . . , α Ng ) ∈ N Ng } be the family of multivariate orthonormal polynomials with respect to a given arbitrary probability measure p Ξ (ξ) dξ on R Ng such that, for all α and β in N Ng , R Ng Ψ α (ξ) Ψ β (ξ) p Ξ (ξ) dξ = E{Ψ α (Ξ) Ψ β (Ξ)} = δ αβ . For α = 0, Ψ 0 (η) = 1 is the constant normalized polynomial. The second-order random vector H ν can then be expanded in polynomial chaos Ψ α as H ν = lim Ng →n,N d →+∞ H ν,N d ,Ng with convergence in the vector space of all the R n -valued second-order random variables, in which H ν,N d ,Ng = K k=0 h (k) Ψ α (k) (Ξ) where h (0) , . . . , h (K) are the coefficients in R n . The statistical inverse problem consists in estimating integers N g ≤ n and N d ≥ 1, and the coefficients h (0) , . . . , h (K) in R n , by using the ν "experimental" realizations η exp,1 , . . . , η exp,ν in order that E{ H ν -H ν,N d ,Ng 2 } ≤ ε E{ H ν 2 }. An adapted method such as those proposed in [START_REF] Arnst | Identification of Bayesian posteriors for coefficients of chaos expansion[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Das | Asymptotic sampling distribution for polynomial chaos representation of data : A maximum-entropy and fisher information approach[END_REF][START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF][START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems[END_REF][START_REF] Nouy | Random fields representations for stochastic elliptic boundary value problems and statistical inverse problems[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data[END_REF] can be used for solving this statistical inverse problem for which the numerical effort is directly related to the value of N g ≤ n. If an optimal partition H ν = (Y ν,1 , . . . , Y ν,m ) is constructed with m ≤ n, then the PCE can be identified for each R µ j -valued random vector Y ν,j such that Y ν,j,N (j) d = K j k=0 h (j,k) Ψ α (j,k) (Ξ j ) in which h (j,0) , . . . , h (j,K j ) are the coefficients in R µ j , where Ξ j is a random variable with values in R N (j) g with N (j) g ≤ µ j < n for which its arbitrary probability measure p Ξ j (ξ j ) dξ j on R N (j) g is given, and where

K j = (N (j) g + N (j) d )!/N (j) g ! N (j) d ! -1.
As the random vectors Ξ 1 , . . . , Ξ m are mutually independent, if m (with m ≤ n) is big, then the gain obtained is very big, because m j=1 K j ≪ K.

(iii) Brief state of the art concerning the methodologies for testing independence.

(1) A popular method for testing the statistical independence of the N components of a random vector from a given set of ν realizations is the use of the frequency distribution [START_REF] Fleiss | Statistical Methods for Rates and Proportions[END_REF] coupled with the use of the Pearson chi-squared (χ 2 ) test [START_REF] Pearson | On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling[END_REF][START_REF] Greenwood | A Guide to Chi-Squared Testing[END_REF]. For the high dimension (N big) and for a relatively small value of ν, such a an approach does not give sufficiently good results. In addition, if this type of method allows for testing the independence, we always need an additional fast algorithm for constructing the optimal partition that we are looking for.

(2) The independent component analysis (ICA) [START_REF] Herault | Space or time adaptive signal processing by neural network models[END_REF][START_REF] Jutten | Blind separation of sources, part I: An adaptive algorithm based on neuromimetic architecture[END_REF][START_REF] Comon | Blind separation of sources, part II: Problems statement[END_REF][START_REF] Comon | Independent component analysis -a new concept?[END_REF][START_REF] Ärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF][START_REF] Comon | Handbook of blind source separation[END_REF][START_REF] Hyv | Independent component analysis: algorithms and applications[END_REF][START_REF] Lee | A unifying information-theoretic framework for independent component analysis[END_REF][START_REF] Bolosco | Independent component analysis based on nonparametric density estimation[END_REF], which is also referred to as the blind source separation, is an efficient method that consists in extracting independent source signals from a linear mixture of mutually statistically independent signals, which is used for source-separation problems, and which is massively used in signal analysis for analyzing, for instance, financial time series or damage in materials, for image processing, in particular for analyzing neuroimage data such as electroencephalogram (EEG) images, neuroimaging of the brain in computational neuroscience, data compression of spectroscopic data sets, etc. The fundamental hypothesis that is introduced in such an ICA methodology is that the observed vector-valued signal is a linear transformation of statistically independent real-valued signals (that is to say, is a linear transformation of a vector-valued signal whose components are mutually independent) and the objective of the ICA algorithms is to identify the best linear operator. The ICA method can be formulated as follows. Let T = {t 1 , . . . , t N T } be a finite set (corresponding to a time index or any other index of interest). For all t in T , let z exp,1 (t), . . . , z exp,ν (t) be ν observations (the realizations) of a time series {Z(t) = (Z 1 (t), . . . , Z N (t)), t ∈ T } with values in R N . In the ICA methodology, it is assumed that, for all t in T , the random vector Z(t) can be modeled as a linear combination of hidden mutually independent random variables Y 1 (t), . . . , Y m (t), with some unknown coefficients

[c] ij , that is to say, Z i (t) = m j=1 [c] ij Y j (t), or Z(t) = [c] Y(t)
, for which matrix [c] and the time series {Y 1 (t), t ∈ T }, . . . , {Y m (t), t ∈ T }, which are assumed to be mutually independent, must be estimated using only {z exp,1 (t), t ∈ T }, . . . , {z exp,ν (t), t ∈ T }.

(3) In this paper, the method proposed and its objective are different from the ICA. The common building block with the existing methods developed for ICA is the use of the mutual information for measuring the statistical independence. All the other parts of the methodology presented in this paper are different and have been constructed in order to obtain a robust construction of an optimal partition of a non-Gaussian random vector in high dimension, which is defined by a relatively small number of realizations.

(iv) Organization of the paper and the methodology proposed. The problem related to the construction of an optimal partition of a non-Gaussian random vector in terms of several mutually independent random vectors is a difficult problem for a random vector in high dimension. The given data set is made up of ν realizations of a non-Gaussian second-order random vector X with values in R N , for which its probability distribution is unknown.

In Section 2, a principal component analysis of X is performed and the ν realizations of the random vector H ν with values in R n with n ≤ N is carried out. Then, the random vector H ν is defined by its probability distribution p H ν (η) dη in which the pdf p H ν is constructed by the kernel density estimation method using the ν realizations of H ν . Section 3 is devoted to the proposed algorithm for constructing an optimal partition of random vector H ν in terms of mutually independent random vectors Y ν,1 , . . . , Y ν,m . For testing the mutual independence of the random vectors Y ν,1 , . . . , Y ν,m of a given partition of random vector H ν , a numerical criterion is defined in coherence with the statistical estimates made with the ν realizations. We have chosen to construct this numerical criterion on the base of the mutual information that is expressed as a function of the Shannon entropy (instead of using another approach such as the one based on the use of the conditional pdf). This choice is motivated by the two following reasons: (i) The entropy of a pdf is expressed as the mathematical expectation of the logarithm of this pdf.

Consequently, for a pdf defined on a high-dimension set (n big), the presence of the logarithm allows for solving the numerical problems induced by the constant of normalization of the pdf.

(ii) As the entropy is a mathematical expectation, this mathematical expectation can be estimated by the Monte Carlo method for which the convergence rate is independent of the dimension thanks to the law of large numbers, a property that is absolutely needed for the high dimensions.

(iii) The criterion that will be constructed by using the Monte Carlo estimation of the entropy will be coherent with the maximum likelihood method, which is often use for solving statistical inverse problems (such as the identification of the coefficients of the PCE of the random vector from a set of realizations).

(iv) We then formulate the problem related to the identification of the optimal partition as an optimization problem. However, the exploration of all the possible partitions that would be necessary for solving this optimization problem is tricky in high dimension. This optimization problem is then replaced by an equivalent one that can be solved with a fast algorithm from the graph theory (detailed in Section 3.6). Section 4 deals with four numerical experiments, which allow for obtaining numerical validations of the approach proposed.

Construction of a non-Gaussian reduced-order statistical model from a set of realizations.

2.1. Data description and estimates of the second-order moments. Let X = (X 1 , . . . , X N ) be a R N -valued non-Gaussian second-order random vector defined on a probability space (Θ, T , P), for which its probability distribution is represented by an unknown pdf x → p X (x) with respect to the Lebesgue measure dx on R N . It is assumed that ν (with ν ≫ 1) independent realizations x exp,1 , . . . , x exp,ν of X are given (coming from experimental data or from numerical simulations). Let m ν X and [ C ν X ] be the empirical estimates of the mean vector m X = E{X} and of the covariance matrix

[C X ] = E{(X -m X ) (X -m X ) T }, such that m ν X = 1 ν ν ℓ=1 x exp,ℓ , [ C ν X ] = 1 ν -1 ν ℓ=1 (x exp,ℓ -m ν X ) (x exp,ℓ -m ν X ) T . (2.1)
We introduce the non-Gaussian second-order random vector X ν with values in R N , defined on (Θ, T , P), for which the ν independent realizations are x exp,1 , . . . , x exp,ν ,

X ν (θ ℓ ) = x exp,ℓ ∈ R N , θ ℓ ∈ Θ , ℓ = 1, . . . , ν , (2.2)
and for which the mean vector and the covariance matrix are exactly m ν X and [ C ν X ],

E{X ν } = m ν X , E{(X ν -E{X ν }) (X ν -E{X ν }) T } = [ C ν X ] . (2.3)

2.2.

Reduced-order statistical model X (n,ν) of X ν . Let n be the reduced-order dimension such that n ≤ N , which is defined (see hereinafter) by analyzing the convergence of the principal component analysis of random vector X ν . The eigenvalues λ 1 ≥ . . . ≥ λ N ≥ 0 and the associated orthonormal eigenvectors ϕ 1 , . . . , ϕ N , such that (ϕ i ) T ϕ j = δ ij , are such that [ C ν X ] ϕ i = λ i ϕ i . The principal component analysis allows for constructing a reduced-order statistical model X (n,ν) of X ν such that

X (n,ν) = m ν X + n i=1 λ i H ν i ϕ i , (2.4)
in which

H ν i = (X ν -m ν X ) T ϕ i / λ i , E{H ν i } = 0 , E{H ν i H ν j } = δ ij . (2.5)
It should be noted that the second-order random variables H ν 1 , . . . , H ν n are non-Gaussian, centered, and uncorrelated, but are statistically dependent. Let (n, ν) → err(n, ν) be the error function defined on {1, . . . , N } × N * such that

err(n, ν) = 1 - n i=1 λ i tr[ C ν X ] , (2.6)
in which λ i depends on ν. Note that n → err(n, ν) depends on ν, but that err(N, ν) = 0 for all ν. For given ε > 0 (sufficiently small with respect to 1), it is assumed that

∃ ν p ∈ N * , ∃ n ∈ {1, . . . , N } : ∀ ν ≥ ν p , err(n, ν) ≤ ε , (2.7)
in which n and ε are independent of ν (for ν ≥ ν p ). In the following, it is assumed that n is fixed such that (2.7) is verified. Consequently, since X (N,ν) = X ν , it can be deduced that

∀ ν ≥ ν p , E{ X ν -X (n,ν) 2 } ≤ ε tr[ C ν X ] . (2.8)
The left-hand side of (2.8) represents the square of the norm of random vector X ν -X (n,ν) in L 2 (Θ, R N ) and consequently, allows for measuring the distance between X ν and X (n,ν) in L 2 (Θ, R N ).

Note that the dimension reduction constructed by using a principal component analysis is not systematically mandatory but is either required for constructing a reduced model such as for a random field or is recommended in order to avoid some possible numerical difficulties that could occur during the computation of the mutual information for which the Gaussian kernel density estimation method is used and therefore, introduces the computation of exponentials.

2.3.

Probabilistic construction of the random vector H ν . In this section, we give the probabilistic construction of the non-Gaussian random vector H ν whose components are the coordinates H ν 1 , . . . , H ν n introduced in the reduced-order statistical model X (n,ν) of X ν , defined by (2.4).

2.3.1.

Independent realizations and second-order moments of the random vector H ν . Let H ν = (H ν 1 , . . . , H ν n ) be the R n -valued non-Gaussian second-order random variable defined on (Θ, T , P). The ν independent realizations {η exp,ℓ = (η exp,ℓ 1 , . . . , η exp,ℓ n ) ∈ R n , ℓ = 1, . . . , ν} of the R n -valued random variable H ν are computed by

η exp,ℓ i = 1 √ λ i (x exp,ℓ -m ν X ) T ϕ i , ℓ = 1, . . . , ν , i = 1, . . . , n , (2.9)
and depend on ν, and for all fixed ν,

m ν H = 1 ν ν ℓ=1 η exp,ℓ = 0 , [ R ν H ] = 1 ν -1 ν ℓ=1 η exp,ℓ (η exp,ℓ ) T = [ I n ] .
(2.10) Equations (2.5) and (2.10) show that

E{H ν } = m ν H = 0 , E{H ν (H ν ) T } = [ R ν H ] = [I n ] .
(2.11)

2.3.2.

Probability distribution of the non-Gaussian random vector H ν . At this step of the construction of random vector H ν , its ν realizations {η exp,ℓ , ℓ = 1, . . . , ν} (see (2.9)) and its second-order moments (see (2.11)) are defined. For completing the probabilistic construction of H ν , we define its probability distribution by the pdf η → p H ν (η) on R n that is constructed by the Gaussian kernel density estimation method on the basis of the knowledge of the ν independent realizations η exp,1 , . . . , η exp,ν . The modification proposed in [START_REF] Soize | Polynomial chaos expansion of a multimodal random vector[END_REF] of the classical Gaussian kernel density estimation method [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis[END_REF] is used so that, for all ν, (2.11) is preserved. The positive valued function p H ν on R n is then written as

p H ν (η) = c n q ν (η) , ∀ η ∈ R n , (2.12)
in which the positive constant c n and the positive-valued function η → q ν (η) on R n are such that

c n = 1 ( √ 2π s n ) n , q ν (η) = 1 ν ν ℓ=1 exp{- 1 2 s 2 n s n s n η exp,ℓ -η 2 } , (2.13)
and where the positive parameters s n and s n are defined by

s n = 4 ν(2 + n) 1/(n+4) , s n = s n s 2 n + ν-1 ν . (2.14)
Parameter s n is the usual multidimensional optimal Silverman bandwidth (taking into account that the standard deviation of each component of H ν is unity) and parameter s n has been introduced in order that the second equation of (2.11) be verified. It should be noted that, for n fixed, parameters s n and s n go to 0 + , and s n /s n goes to 1 -, when ν goes to +∞. Using (2.12) to (2.14), it can easily be verified that

E{H ν } = R n η p H ν (η) dη = s n s n m ν H = 0 , (2.15
)

E{H ν (H ν ) T } = R n η η T p H ν (η) dη = s 2 n [I n ] + s n s n 2 (ν -1) ν [ R ν H ] = [ I n ] . (2.16)
Remark. By construction of p H ν , the sequence {p H ν } ν has a limit, denoted by p H = lim ν→+∞ p H ν , in the space of all the integrable positive-valued functions on R n . Consequently, the sequence {H ν } ν of R n -valued random variables has a limit, denoted by H = lim ν→+∞ H ν , in probability distribution. The probability distribution of the R n -valued random variable H is defined by the pdf p H .

2.3.3.

Marginal distribution related to a subset of the components of random vector H ν . Let j 1 < j 2 < . . . < j µ j with µ j < n be a subset of {1, 2, . . . , n}. Let Y ν,j = (H j 1 , H j 2 , . . . , H jµ j ) be the random vector with values in R µ j and let be η j = (η j 1 , . . . , η j µ j ) ∈ R µ j . From (2.12) and (2.13), it can easily be deduced that the pdf p Y ν,j of Y ν,j with respect to dη j , which is such that p Y ν,j (η j ) = p H j 1 ,H j 2 ,...,H jµ j (η j 1 , . . . , η j µ j ), can be written as

p Y ν,j (η j ) = c µ j q ν j (η j ) , ∀ η j ∈ R µ j , (2.17)
in which the positive constant c µ j and the positive-valued function η j → q ν j (η j ) on R µ j are such that

c µ j = 1 ( √ 2π s n ) µ j , q ν j (η j ) = 1 ν ν ℓ=1 exp{- 1 2 s 2 n µ j k=1 ( s n s n η exp,ℓ j k -η j k ) 2 } . (2.18)
Note that, in (2.18), s n and s n must be used and not s µ j and s µ j .

Remark concerning convergence aspects.

It can be seen that, for ν → +∞, the sequence of random vectors {X ν } ν converges in probability distribution to random vector X for which the probability distribution P X (dx) = p X (x) dx on R N is the limit (in the space of the probability measures on R N ) of the sequence of the probability distributions {P X ν (dx) = p X ν (x) dx} ν in which p X ν is the pdf of the random vector X ν = X (N,ν) given by (2.4) with n = N , and where the pdf p H ν of random vector H ν is defined by (2.12). It should be noted that this result holds because it is assumed that X is a second-order random vector and that its probability distribution P X (dx) admits a density p X with respect to dx.

3. An algorithm for constructing an optimal partition of random vector H ν in terms of mutually independent random vectors Y ν,1 , . . . , Y ν,m . For n fixed such that n ≤ N and for any ν fixed such that ν ≥ ν p , as H ν = (H ν 1 , . . . , H ν n ) is a normalized random vector (centered and with a covariance matrix equal to [I n ], see (2.11)), if H ν was Gaussian, then the components H ν 1 , . . . , H ν n would be independent. As H ν is assumed to be a non-Gaussian random vector, although H ν is normalized, its components H ν 1 , . . . , H ν n are, a priori, mutually dependent (it could also be non-Gaussian and mutually independent). In this section that is central for this paper, the question solved is related to the development of an algorithm for constructing an optimal partition of H ν , which consists in finding the largest value m max ≥ 1 of the number m of random vectors Y ν,1 , . . . , Y ν,m made up of the components H ν 1 , . . . , H ν n , such that H ν can be written as

H ν = (Y ν,1 , . . . , Y ν,m
) in which the random vectors Y ν,1 , . . . , Y ν,m are mutually independent, but for which the components of each random vector Y ν,j are, a priori, mutually dependent. If the largest value m max of m is such that

• m max = 1, then there is only one random vector, Y ν,1 = H ν , and all the components of H ν are mutually dependent and cannot be separated into several random vectors; • m max = n, then each random vector of the partition is made up of one component of H ν , which means that, Y ν,j = H ν j and all the components H ν 1 , . . . , H ν n of H ν are mutually independent.

As we have explained in Section 1, for a random vector H ν in high dimension (big value of n) defined by a set of realizations, this problem is difficult. An efficient algorithm based on the use of information theory and of graph theory is proposed hereinafter.

Defining a partition of H

ν in terms of mutually independent random vectors. Let m be an integer such that 1 ≤ m ≤ n .

A partition of the components H ν 1 , . . . , H ν n of H ν , denoted by P ν (m; µ 1 , . . . , µ m ), consists in writing H ν in terms of m random vectors Y ν,1 , . . . , Y ν,m such that

H ν = (Y ν,1 , . . . , Y ν,m ) , (3.1)
in which, for all j in {1, . . . , m}, µ j is an integer such that 1 ≤ µ j ≤ n and Y ν,j = (Y ν,j 1 , . . . , Y ν,j µ j ) is a R µ j -valued random vector, which is written as

Y ν,j = (H ν r j 1 , . . . , H ν r j µ j ) , 1 ≤ r j 1 < r j 2 < . . . < r j µ j ≤ n , (3.2)
in which the equality r j µ j = n can exist only if m = 1 (and thus

µ 1 = n that yields Y ν,1 = H ν ). For fixed ν, the integers µ 1 ≥ 1, . . . , µ m ≥ 1 are such that ∪ m j=1 {r j 1 , . . . , r j µ j } = {1, . . . , n} , ∩ m j=1 {r j 1 , . . . , r j µ j } = {∅} , µ 1 + . . . + µ m = n . (3.3)
It should be noted that (3.2) implies that there is no possible repetition of indices r j 1 , . . . , r j µ j and that no permutation is considered (any permutation would yield the same probability distribution for

Y ν,j ). If the partition P ν (m; µ 1 , . . . , µ m ) is made up of m mutually independent random vectors Y ν,1 , . . . , Y ν,m , then the joint pdf (η 1 , . . . , η m ) → p Y ν,1 ,...,Y ν,m (η 1 , . . . , η m ) defined on R µ 1 × . . . × R µm is such that, for all η = (η 1 , . . . , η m ) in R n = R µ 1 × . . . × R µm , p H ν (η) = p Y ν,1 ,...,Y ν,m (η 1 , . . . , η m ) = p Y ν,1 (η 1 ) × . . . × p Y ν,m (η m ) , (3.4)
where, for all j in {1, . . . , m}, η j → p Y ν,j (η j ) is the pdf on R µ j of random vector Y ν,j . From (2.12) and (3.1), it can be deduced that the joint pdf p

Y ν,1 ,...,Y ν,m on R µ 1 × . . . × R µm of the random vectors Y ν,1 , . . . , Y ν,m can be written, for all η = (η 1 , . . . , η m ) in R n = R µ 1 × . . . × R µm , as p Y ν,1 ,...,Y ν,m (η 1 , . . . , η m ) = p H ν (η) = c n q ν (η) , (3.5)
in which c n and q ν are defined by (2.13).

Setting the problem for finding an optimal partition P ν

opt in terms of mutually independent random vectors. For n and ν fixed such that (2.7) is satisfied, the problem related to the construction of an optimal partition P ν opt = P ν (m max ; µ opt 1 , . . . , µ opt mmax ) in terms of mutually independent random vectors Y ν,1 , . . . , Y ν,mmax consists in finding the largest number m max of integer m on the set of all the possible partitions P ν (m; µ 1 , . . . , µ m ). As we will explain later, a numerical criterion must be constructed for testing the mutual independence of the random vectors Y ν,1 , . . . , Y ν,m in a context for which, in general, the number ν of the given realizations of H ν is not sufficiently large so that the statistical estimators are well converged (negligible statistical fluctuations of the estimates).

Construction of a numerical criterion for testing a partition in terms of mutually independent random vectors.

In order to characterize the mutual independence of the random vectors Y ν,1 , . . . , Y ν,m of a partition P ν (m; µ 1 , . . . , µ m ) of random vector H ν , we need to introduce a numerical criterion that we propose to construct by using the mutual information [START_REF] Kolmogorov | On the Shannon theory of information transmission in the case of continuous signals[END_REF][START_REF] Cover | Elements of Information Theory[END_REF]), which is defined as the relative entropy (introduced by Kullback and Leibler [START_REF] Leibler | On information and sufficiency[END_REF]) between the joint pdf of all the random vectors Y ν,1 , . . . , Y ν,m (that is to say, the pdf of H ν ) and the product of the pdf of each random vector Y ν,j . This mutual information and the relative entropy can be expressed with the entropy (also called the differential entropy or the Shannon entropy, introduced by Shannon [START_REF] Shannon | A MATHEMATICAL THEORY OF COMMUNICATION[END_REF]) of the pdf of the random vectors. This choice of constructing a numerical criterion based on the use of the entropy (instead of using another approach such as the one based on the use of the conditional pdf) is motivated by the following reasons:

• As it has previously been explained, we recall that the entropy of a pdf is expressed as the mathematical expectation of the logarithm of this pdf. Consequently, for a pdf whose support is a set having a high dimension, the presence of the logarithm allows for avoiding the numerical problems that are induced by the presence of the constant of normalization of the pdf.

• As the entropy is a mathematical expectation, this mathematical expectation can be estimated by the Monte Carlo method for which the convergence rate is independent of the dimension thanks to the law of large numbers (central limit theorem) [START_REF] Givens | Computational Statistics[END_REF][START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF][START_REF] Serfling | Approximation Theorems of Mathematical Statistics[END_REF], property that is absolutely needed for the high dimensions.

• The numerical criterion that will be constructed by using the Monte Carlo estimation of the entropy will be coherent with the maximum likelihood method, which is, for instance (see Section 1), used for the identification of the coefficients of the PCE of random vector H ν by using the set of realizations η exp,1 , . . . , η exp,ν [START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data[END_REF]. Nevertheless, a numerical criterion that is based on the use of the mutual information cannot be constructed as a direct application of the mutual information, and some additional specific ingredients must be introduced for obtaining a numerical criterion that is robust with respect to the statistical fluctuations induced by a nonperfect convergence of the statistical estimators used, because, for many applications, the number ν of the available independent realizations of the random vector H ν is not sufficiently big.

Entropy of the pdfs related to a partition.

For all j in {1, . . . , m}, the entropy of the pdf p Y ν,j of the R µ j -valued random variable Y ν,j is S(p Y ν,j ) ∈ R that is defined (using the convention 0 log(0) = 0) by

S(p Y ν,j ) = -E{log(p Y ν,j (Y ν,j ))} = - R µ j p Y ν,j (η j ) log(p Y ν,j (η j )) dη j , (3.6)
in which pdf p Y ν,j is defined by (2.17). Similarly, the entropy S(p

Y ν,1 ,...,Y ν,m ) ∈ R of the pdf p Y ν,1 ,...,Y ν,m is defined by S(p Y ν,1 ,...,Y ν,m ) = -E{log(p Y ν,1 ,...,Y ν,m (Y ν,1 , . . . , Y ν,m ))} , (3.7)
which can be rewritten (using (3.5)) as

S(p Y ν,1 ,...,Y ν,m ) = S(p H ν ) = -E{log(p H ν (H ν ))} = - R n p H ν (η) log(p H ν (η)) dη . (3.8)

Mutual information of the random vectors related to a partition. The mutual information

i(Y ν,1 , . . . , Y ν,m ) between the random vectors Y ν,1 , . . . , Y ν,m is defined by i(Y ν,1 , . . . , Y ν,m ) = -E log( p Y ν,1 (Y ν,1 ) × . . . × p Y ν,m (Y ν,m ) p Y ν,1 ,...,Y ν,m (Y ν,1 , . . . , Y ν,m ) ) , (3.9
Y ν,1 ⊗ . . . ⊗ p Y ν,m , d(p Y ν,1 ,...,Y ν,m : p Y ν,1 ⊗ . . . ⊗ p Y ν,m ) = i(Y ν,1 , . . . , Y ν,m ) . (3.10)
In the following, we will use the mutual information that is equivalent to the Kullback-Leibler divergence. From (3.6) and (3.7), it can be deduced that the mutual information defined by (3.9) can be written as

i(Y ν,1 , . . . , Y ν,m ) = S(p Y ν,1 ) + . . . + S(p Y ν,m ) -S(p Y ν,1 ,...,Y ν,m ) . (3.11)

Inequalities verified by the mutual information and by the entropy.

We have the classical important following inequalities that allow for testing the mutual independence of the random vectors defined by a partition P ν (m; µ 1 , . . . , µ m ): (i) Using the Jensen inequality, it can be proved [START_REF] Cover | Elements of Information Theory[END_REF] that the mutual information is non negative and can be equal to +∞, Taking into account the properties of the mutual information, for testing the mutual independence of random vectors Y ν,1 , . . . , Y ν,m of a partition P ν (m; µ 1 , . . . , µ m ), we choose the mutual information

0 ≤ i(Y ν,1 , . . . , Y ν,m ) ≤ +∞ . (3.12)
0 ≤ i(Y ν,1 , . . . , Y ν,m ) ≤ +∞ as the theoretical criterion that is such that i(Y ν,1 , . . . , Y ν,m ) = 0 if Y ν,1 , . . . , Y ν,m are mutually independent. Equations (2.13) and (2.18) yield c n = c µ 1 × . . . × c µm .
Consequently, (2.12) and (2.17) show that the mutual information i(Y ν,1 , . . . , Y ν,m ), which is defined by (3.11) and deduced from (3.9), is independent of the constants of normalization c n , c µ 1 , . . . , c µm . However, such a theoretical criterion cannot directly be used (and will be replaced by a numerical criterion derived from the theoretical criterion) for the following reasons:

• The number of realizations ν is fixed and can be relatively small (this means that the asymptotic value ν → +∞ is not reached at all for obtaining a good convergence of the statistical estimators). Consequently, there are statistical fluctuations of the theoretical criterion as a function of ν, which exclude the possibility for obtaining exactly the value 0 (by superior values). ) cannot exactly be calculated. As we have explained before, for the high dimensions, the Monte Carlo method must be used and, consequently, only a numerical approximation can be obtained. As p H ν is known, a generator of realizations could be used for generating a large number of independent realizations in order to obtain a very good approximation of these entropies. Such a generator could be constructed by using a Markov Chain Monte Carlo (MCMC) algorithm such as the popular Metropolis-Hastings algorithm [START_REF] Givens | Computational Statistics[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF][START_REF] Metropolis | The Monte Carlo method[END_REF], the Gibbs sampling [START_REF] Geman | Stochastic relaxation, Gibbs distribution and the Bayesian distribution of images[END_REF], or other algorithms. Nevertheless, in a following section, for estimating these entropies, we will propose to use only the given set of realizations η exp,1 , . . . , η exp,ν without adding some additional realizations computed with an adapted generator. The reason why is we want to use exactly the same data for testing the independence and for constructing the estimators related to the identification of the representation of H ν under consideration in order to preserve the same level of statistical fluctuations in all the estimators used.

Remarks concerning the choice of the theoretical criterion. The theoretical criterion consists in testing the mutual information i(Y ν,1 , . . . , Y ν,m ) with respect to zero (by superior values) because i(Y ν,1 , . . . , Y ν,m ) ≥ 0, and the mutual independence is obtained for 0. As this type of test (with respect to 0) has no numerical meaning, this theoretical criterion must be replaced by a numerical criterion that will be derived, and which will be presented in Section 3.3.5. However, it could be interesting to examine what could be other choices for defining a theoretical criterion for mutual independence. Hereinafter, we examine two possibilities among others: (i) Considering (3.11), (3.13), and (3.14), a first alternative choice would consist in introducing the quantity S(p H ν )/ m j=1 S(p Y ν,j ) ≤ 1 that is equal to 1 when Y ν,1 , . . . , Y ν,m are mutually independent, but which requires us to introduce the hypothesis:

m j=1 S(p Y ν,j ) > 0.
Unfortunately, such a hypothesis is not satisfied in the general case. (ii) A second choice could be introduced on the base of the following result. Let G ν be the Gaussian second-order centered R n -valued random vector for which its covariance matrix is [I n ]. Consequently, its components are mutually independent and its entropy is S(p G ν ) = n (1 + log(2π))/2 > 1. On the other hand, it is known (see [START_REF] Cover | Elements of Information Theory[END_REF]) that for any non-Gaussian centered random vector H ν with covariance matrix equal to

[I n ], we have S(p H ν ) ≤ S(p G ν ). It should be noted that the equality S(p H ν ) = S(p G ν ) is verified only if H ν is Gaussian.
The components of H ν can be independent and not Gaussian, and in such a case, we have S(p H ν ) < S(p G ν ). Thus, for all j = 1, . . . , m, we have S(p Y ν,j ) ≤ S(p G ν,j ) and it can then be deduced that 0

≤ m j=1 S(p Y ν,j ) -S(p H ν ) ≤ m j=1 S(p G ν,j ) -S(p H ν ) = S(p G ν ) -S(p H ν ). For the non-Gaussian case, S(p G ν ) -S(p H ν ) > 0 and, consequently, the theoretical criterion 0 ≤ { m j=1 S(p Y ν,j ) -S(p H ν )}/{S(p G ν ) -S(p H ν )} ≤ 1 could be constructed.
Unfortunately, such a criterion would be ambiguous because 0 (by superior values) corresponds to the independence of non-Gaussian vectors, but for a small non-Gaussian perturbation of a Gaussian vector, that is to say if

H ν → G ν in probability distribution, then S(p G ν ) -S(p H ν ) → 0
(by superior values) and, consequently, an indeterminate form of the criterion is, a priori, obtained.

Numerical criterion for testing the mutual independence of the random vectors of a partition.

As we have explained in Section 3.3.4, the theoretical criterion cannot directly be used and a numerical criterion must be derived. Let us introduce the positive-valued random variable Z ν such that

Z ν = p Y ν,1 (Y ν,1 ) × . . . × p Y ν,m (Y ν,m ) p Y ν,1 ,...,Y ν,m (Y ν,1 , . . . , Y ν,m ) . (3.16)
The theoretical criterion (mutual information) defined by (3.9) can be rewritten as

i(Y ν,1 , . . . , Y ν,m ) = -E{log Z ν } . (3.17)
(i) Defining an approximate criterion for testing the mutual independence of the random vectors of a partition. An approximation of the theoretical criterion defined by (3.9), and rewritten as (3.17), is constructed by introducing the usual estimator of the mathematical expectation, but, as we have previously explained, by using only the ν independent realizations η exp,1 , . . . , η exp,ν of H ν . As η = (η 1 , . . . , η m ) ∈ R n = R µ 1 ×. . .×R µm is the deterministic vector associated with random vector H ν = (Y ν,1 , . . . , Y ν,m ), for ℓ = 1, . . . , ν, the realization η exp,ℓ ∈ R n is rewritten as η exp,ℓ = (η 1,exp,ℓ , . . . , η m,exp,ℓ ) with η j,exp,ℓ = (η j,exp,ℓ 1 , . . . , η j,exp,ℓ

µ j ) = (η exp,ℓ r j 1 , . . . , η exp,ℓ r j µ j
) ∈ R µ j . Let Z exp,1 , . . . , Z exp,ν be ν independent copies of random variable Z ν . The approximate criterion for testing the mutual independence of random vectors Y ν,1 , . . . , Y ν,m is defined as the real-valued random variable

I ν (Y ν,1 , . . . , Y ν,m ) = - 1 ν ν ℓ=1 log Z exp,ℓ . (3.18)
If the random vectors Y ν,1 , . . . , Y ν,m are mutually independent, then I ν (Y ν,1 , . . . , Y ν,m ) = 0 almost surely.

(ii) Computing a realization i ν (Y ν,1 , . . . , Y ν,m ) of the random variable I ν (Y ν,1 , . . . , Y ν,m ). The realization of I ν (Y ν,1 , . . . , Y ν,m ) associated with the ν independent realizations η exp,1 , . . . , η exp,ν of H ν is denoted by i ν (Y ν,1 , . . . , Y ν,m ), which can be written (see (3.6), (3.8), (3.11), and (3.16) to (3.18)) as

i ν (Y ν,1 , . . . , Y ν,m ) = s ν 1 + . . . + s ν m -s ν , (3.19)
in which the real numbers {s ν j , j = 1, . . . , m} and s ν can be computed with the formulas,

s ν j = - 1 ν ν ℓ=1 log(p Y ν,j (η j,exp,ℓ )) , j = 1, . . . , m , (3.20
)

s ν = - 1 ν ν ℓ=1 log(p H ν (η exp,ℓ )) . (3.21)
If the random vectors Y ν,1 , . . . , Y ν,m are mutually independent, then i ν (Y ν,1 , . . . , Y ν,m ) = 0.

(iii) Connection with the likelihood function. From (3.20) and (3.21), s ν j and s ν can be rewritten with the log likelihood functions L ν j = ν ℓ=1 log(p Y ν,j (η j,exp,ℓ )) and L ν = ν ℓ=1 log(p H ν (η exp,ℓ )) as s ν j = -L ν j /ν and s ν = -L ν /ν. The realization of the approximate criterion can then be rewritten as

i ν (Y ν,1 , . . . , Y ν,m ) = -(L ν 1 + . . . + L ν m -L ν )/ν.
(iv) Random lower bound of the approximate criterion I ν (Y ν,1 , . . . , Y ν,m ). As the independent random variables Z exp,1 , . . . , Z exp,ν are positive almost surely, and as the arithmetic mean (Σ ν ℓ=1 Z exp,ℓ )/ν is greater than or equal to the geometric mean (Π ν ℓ=1 Z exp,ℓ ) 1/ν , it can be deduced that

log{(Σ ν ℓ=1 Z exp,ℓ )/ν} ≥ log{(Π ν ℓ=1 Z exp,ℓ ) 1/ν } = 1 ν ν ℓ=1 log(Z exp,ℓ ) a.s . (3.22)
Consequently, (3.18) and (3.22) yield the following inequality of real-valued random variables,

B ν ≤ I ν (Y ν,1 , . . . , Y ν,m ) a.s , (3.23)
in which the real-valued random variable B ν is written as

B ν = -log( 1 ν ν ℓ=1 Z exp,ℓ ) . (3.24)
If the random vectors Y ν,1 , . . . , Y ν,m are mutually independent, then B ν = 0 almost surely. Considering the ν independent realizations η exp,1 , . . . , η exp,ν of H ν , (3.16) and (3.24) yield the deterministic inequality

0 ≤ b ν ≤ i ν (Y ν,1 , . . . , Y ν,m ) , (3.25)
in which the real number b ν is given by

b ν = -log 1 ν ν ℓ=1 p Y ν,1 (η 1,exp,ℓ ) × . . . × p Y ν,m (η m,exp,ℓ ) p H ν (η exp,ℓ ) . (3.26)
If the random vectors Y ν,1 , . . . , Y ν,m are mutually independent, then b ν = 0. It should be noted that ν p is assumed to be such that, for all ν ≥ ν p , we have effectively 0 ≤ b ν , which implies that for all ν ≥ ν p , i ν (Y ν,1 , . . . , Y ν,m ) ≥ 0. Unfortunately, as ν is generally not sufficiently large for getting good convergence of the estimator of the pdfs in the right-hand side of (3.26), the lower bound b ν cannot directly be used in the algorithm. An adapted numerical criterion is introduced in paragraph (vi) hereafter.

(v) Convergence analysis. Taking into account (2.7), n is fixed and is independent of ν. From the remark given in Section 2.3.2, the sequence {H ν } ν of random vectors tends in probability distribution to the random vector H when ν goes to +∞. The probability distribution of H is defined by the pdf p H . Similarly to (3.1), the limit H is written as H = (Y 1 , . . . , Y m ) in which for all j in {1, . . . , m}, Y j = (Y j 1 , . . . , Y j µ j ) is an R µ j -valued random vector, where the integers m and µ 1 , . . . , µ m are fixed, independent of ν, and verify (3.3) (see Section 3.1). Consequently, for ν → +∞, i(Y ν,1 , . . . , Y ν,m ) → i(Y 1 , . . . , Y m ). It is assumed that pdf p H (that is unknown) is such that there exists ν 0 such that, for all ν > ν 0 , E{(log Z ν ) 2 } ≤ c Z < +∞ in which the positive constant c Z is independent of ν. For all ν fixed, it can easily be proved that m

I ν = E{I ν (Y ν,1 , . . . , Y ν,m )} = i(Y ν,1 , . . . , Y ν,m ) and σ 2 I ν = E{(I ν (Y ν,1 , . . . , Y ν,m ) -m I ν ) 2 } = 1 ν σ 2 log Z ν with σ 2 log Z ν = E{(log Z ν ) 2 } -i(Y ν,1 , . . . , Y ν,m ) 2 . Consequently, σ 2 I ν → 0 as ν → 0.
For all ε > 0, the use of the Chebyshev inequality yields

P{|I ν (Y ν,1 , . . . , Y ν,m ) -i(Y ν,1 , . . . , Y ν,m )| ≥ ε} ≤ σ 2
I ν /ε 2 that allows for obtaining the following convergence property in probability, It should be noted that the mutual information of the random vectors Y 1 , . . . , Y m is such that

lim ν→+∞ I ν (Y ν,1 , . . . , Y ν,m ) = i(Y 1 , . . . , Y m ) . ( 3
0 ≤ i(Y 1 , . . . , Y m ) , (3.29)
and, if the random vectors Y 1 , . . . , Y m are independent, then i(Y 1 , . . . , Y m ) = 0.

(vi) Defining a numerical criterion for testing the mutual independence of the random vectors of a partition P ν (m; µ 1 , . . . , µ m ). For fixed ν, we have proved that b ν ≤ i ν (Y ν,1 , . . . , Y ν,m ) and, if the random vectors Y ν,1 , . . . , Y ν,m are mutually independent, then

0 = b ν = i ν (Y ν,1 , . . . , Y ν,m ).
Unfortunately, often, ν is not sufficiently large for getting good convergence of the estimator of the pdfs, and consequently, the mutual information i ν (Y ν,1 , . . . , Y ν,m ) cannot be compared to the value 0 for testing the mutual independence. We thus propose to construct a numerical criterion derived from the approximate criterion in order to test the mutual independence, as explained hereinafter.

Let G = (G 1 , . . . , G n ) be the Gaussian second-order centered R n -valued random vector for which its covariance matrix is [I n ]. It is assumed that G is statistically independent of random vector H ν . Therefore the components G 1 , . . . , G n are mutually independent real-valued random variables. Let {g exp,ℓ , ℓ = 1, . . . , ν} be ν independent realizations of G such that g exp,ℓ = (g exp,ℓ 1 , . . . , g exp,ℓ n ) ∈ R n . • By applying the partition P ν (m; µ 1 , . . . , µ m ) defined by (3.1) to (3.3) to random vector G, this vector is rewritten as G = (G 1 , . . . , G m ) in which the Gaussian random vector

G j = (G j r j 1 , . . . , G j r j µ j
) has values in R µ j , where the indices 1 ≤ r j 1 < r j 2 < . . . < r j µ j ≤ n are those defined in (3.2) and (3.3), and where the ν independent realizations of G j are {g j,exp,ℓ , ℓ = 1, . . . , ν} such that g j,exp,ℓ = (g j,exp,ℓ 1 , . . . , g j,exp,ℓ µ j ) = (g exp,ℓ r j 1 , . . . , g exp,ℓ r j µ j

) ∈ R µ j .

• The mutual information i ν (G 1 , . . . , G m ) of random vectors G 1 , . . . , G m is calculated by using (3.19) to (3.21) after replacing random vectors Y ν,1 , . . . , Y ν,m by the mutually independent Gaussian random vectors G 1 , . . . , G m . However, although the values of the pdf, p G (g exp,ℓ ) and p G j (g j,exp,ℓ )), could be calculated exactly, these quantities are approximated by using the Gaussian kernel density estimation method as explained in Sections 2.3.2 and 2.3.3 (in order to define a numerical criterion that is coherent from a numerical point of view).

• The Gaussian random vectors G 1 , . . . , G m are mutually independent. For ν → +∞, b ν Gauss = i ν (G 1 , . . . , G m ) goes to 0. Nevertheless, as ν is finite (and sometimes small enough), the numerical value i ν (G 1 , . . . , G m ) will not be zero, but will be such that 0 < b ν Gauss ≤ i ν (G 1 , . . . , G m ).

Consequently, for ν fixed, the numerical value i ν (G 1 , . . . , G m ) is chosen as the lower bound of the mutual information for testing the mutual independence. • The random vectors Y ν,1 , . . . , Y ν,m of the partition P ν (m; µ 1 , . . . , µ m ) will be considered as mutually independent if

0 ≤ i ν (Y ν,1 , . . . , Y ν,m ) ≤ i ν (G 1 , . . . , G m ) , (3.30)
while they will be considered as mutually dependent if

0 ≤ i ν (G 1 , . . . , G m ) < i ν (Y ν,1 , . . . , Y ν,m ) . (3.31)
• Taking into account (3.25), (3.30), and (3.31), for testing the independence of the random vectors Y ν,1 , . . . , Y ν,m of partition P ν (m; µ 1 , . . . , µ m ), we introduce τ ν (m; µ 1 , . . . , µ m ) as the numerical criterion, which is defined by

τ ν (m; µ 1 , . . . , µ m ) = 1 - i ν (Y ν,1 , . . . , Y ν,m ) i ν (G 1 , . . . , G m ) . (3.32)
This numerical criterion is such that random vectors Y ν,1 , . . . , Y ν,m will be considered as

mutually independent if 0 ≤ τ ν (m; µ 1 , . . . , µ m ) ≤ 1 , mutually dependent if τ ν (m; µ 1 , . . . , µ m ) < 0 .
By convention, τ ν (m; µ 1 , . . . , µ m ) will be taken equal to

1 if i ν (G 1 , . . . , G m ) ≪ 1 with i ν (Y ν,1 , . . . , Y ν,m ) ≤ i ν (G 1 , . . . , G m
). The numerical criterion defined by (3.32) has been constructed for analyzing the high-dimension case (n big) with a number of realizations, which can be relatively small (ν small). Consequently, the estimate i ν (Y ν,1 , . . . , Y ν,m ) of the mutual information performed by using the nonparametric statistics is not sufficiently accurate and there exist statistical fluctuations. In order to increase the robustness of the prevision of the mutual independence, the theoretical lower bound "zero" is replaced by the lower bound b ν Gauss (that would be equal to zero if ν was equal to infinity), and which is estimated with the same number ν of realizations as for i ν (Y ν,1 , . . . , Y ν,m ), in order to construct a robust numerical criterion that takes into account the statistical fluctuations of the lower bound and that allows for testing the mutual independence of the random vectors of a given partition. Finally, it should be noted that the numerical criterion defined by (3.32)

corresponds to a lower bound b ν

Gauss that depends on the partition analyzed and, consequently, the numerical criterion changes for every partition, but, this is the effect that is searched in order to adapt the level of the statistical fluctuations of the numerical criterion to every partition analyzed.

(vii) Comment about the introduction of random variable B ν . It should be noted that the random variable B ν and its estimate b ν have been introduced for justifying the introduction of the lower bound for the Gaussian case, which allows for constructing the numerical criterion defined by (3.32). Consequently, B ν and b ν are introduced for a theoretical understanding of the proposed construction, but do not play any role in the algorithm that will be introduced in Sections 3.5 and 3.6. Remarks about the construction proposed.

Optimization problem for finding the optimal partition of the random vector H

• The optimal partition that is constructed depends on ν. Such a formulation is reasonable because the available information is only made up of the set of the ν realizations of random vector H ν . • The principle that corresponds to the formulation defined by (3.33) consists in finding the partition in terms of mutually independent vectors for which (3.30) is obtained with the largest positive value of the difference i ν (G 1 , . . . , G m )i ν (Y ν,1 , . . . , Y ν,m ), that is to say, with the largest positive value of τ ν (m; µ 1 , . . . , µ m ). Such a criterion yields the largest robustness with respect to the statistical fluctuations induced by the fact that ν can be small enough. • The mathematical analysis of the optimization problem defined by (3.33) seems difficult enough and we have not been able to obtain interesting properties in this way, as against many numerical experiments have confirmed that this choice was efficient with respect to the problem that we had to solve. • For m max ≥ 2, let T = arg max 1≤m≤n {max (µ 1 ,...,µm)∈Cm τ ν (m; µ 1 , . . . , µ m )} be the set of all the integers m with 1 ≤ m ≤ n such that m → max (µ 1 ,...,µm)∈Cm τ ν (m; µ 1 , . . . , µ m ) is maximum with respect to m. A priori, set T is not always reduced to a single point. If m inf max is the number of independent random vectors of the optimal partition corresponding to inf m {T } and m sup max the one corresponding to sup m {T }, then it can be deduced that m inf max ≤ m sup max . Although we are interested in identifying the largest value of m, the most robust solution leads us to choose the largest number of dependent random vectors for the identified optimal partition. Therefore, m max = m inf max = inf m {T } is selected as the maximum number of independent random vectors for this optimal partition. This is the reason why the inferior value of set T has been introduced in (3.33).

Reformulation of the optimization problem for constructing the optimal partition of the random vector H ν in terms of mutually independent random vectors.

The reason why a reformulation of the optimization problem defined by (3.33) is necessary is explained hereinafter, and a new formulation that will be adapted to the development of an algorithm based on graph theory (presented in Section 3.6) is proposed.

Why a reformulation of the optimization problem is necessary.

The optimization problem defined by (3.33) consists in finding the optimal partition P ν opt = P(m max ; µ opt 1 , . . . , µ opt mmax ) on the set of all the possible partitions {P ν (m; µ 1 , . . . , µ m ), 1 ≤ m ≤ n, µ 1 + . . . + µ m = n}, for which µ opt 1 + . . . + µ opt mmax = n and where m max is the largest value of the number m of independent random vectors Y ν,1 , . . . , Y ν,m such that, for j = 1, . . . , m,

Y ν,j = (H ν r j 1 , . . . , H ν r j µ j
) , 1 ≤ r j 1 < r j 2 < . . . < r j µ j ≤ n verifying (3.3) .

(3.34)

For solving (3.33), a natural algorithm would consist in computing τ ν (m; µ 1 , . . . , µ m ) for each partition P ν (m; µ 1 , . . . , µ m ) among all the possible partitions. However, as the number of all the possible partitions is n-1 j=1 C j n with C j n = ((nj + 1) × . . . × n)/j! , it can be seen that, for a high value of n, such a computation would be very tricky. We then propose another algorithm, which is faster than the natural one. This algorithm will be based on the use of the graph theory. The edges of the graph that will be constructed are obtained by analyzing the statistical dependence of the components (H ν r 1 , H ν r 2 ), two by two, for which the indices r 1 and r 2 (with r 1 = r 2 ) are taken in the set of the indices of components H ν 1 , . . . , H ν n of H ν . Consequently, first, we need to adapt the previous notations relative to the mutual information for two scalar random components and, second, we need to reformulate the optimization problem defined by (3.33).

Adapting the notations for the computation of the mutual information of the real-valued random variables H ν r

1 and H ν r 2 . Let r 1 and r 2 be any two integers belonging to {1, . . . , n} with r 1 = r 2 . Hereinafter, the results related to the mutual information and to the numerical criterion presented in Section 3.3 are adapted for analyzing the independence of the real-valued random variables H ν r 1 and H ν r 2 .

Let i ν (H ν r 1 , H ν r 2 ) be the realization of the random mutual information I ν (H ν r 1 , H ν r 2 ), computed by using the ν independent realizations η exp,1 , . . . , η exp,ν of H ν , in which

η exp,ℓ = (η exp,ℓ 1 , . . . , η exp,ℓ n ) ∈ R n . From (3.19) to (3.21), it can be deduced that i ν (H ν r 1 , H ν r 2 ) = s ν,2 r 1 + s ν,2 r 2 -s ν,2 r 1 ,r 2 , (3.35)
in which the real numbers s ν,2 r 1 , s ν,2 r 2 , and s ν,2 r 1 ,r 2 are computed with the formulas,

s ν,2 r = - 1 ν ν ℓ=1 log(p H ν r (η exp,ℓ r
)) , r ∈ {r 1 , r 2 } , (3.36)

s ν,2 r 1 ,r 2 = - 1 ν ν ℓ=1 log(p H ν r 1 ,H ν r 2 (η exp,ℓ r 1 , η exp,ℓ r 2 )) . (3.37)
From Sections 2.3.2 and 2.3.3, it can be seen that the pdf of H ν r is written as

p H ν r (η r ) = c 1 q ν r (η r ) , ∀ η r ∈ R , (3.38)
in which the positive constant c 1 and the positive-valued function η r → q ν r (η r ) on R are such that

c 1 = 1 √ 2π s n , q ν r (η r ) = 1 ν ν ℓ=1 exp{- 1 2 s 2 n ( s n s n η exp,ℓ r -η r ) 2 } . (3.39)
The joint pdf of (H ν r 1 , H ν r 2 ) is written as

p H ν r 1 ,H ν r 2 (η r 1 , η r 2 ) = c 2 q ν r 1 ,r 2 (η r 1 , η r 2 ) , ∀ (η r 1 , η r 2 ) ∈ R 2 , (3.40)
in which the positive constant c 2 and the positive-valued function

(η r 1 , η r 2 ) → q ν r 1 ,r 2 (η r 1 , η r 2 ) on R 2 are such that c 2 = 1 ( √ 2π s n ) 2 , q ν r 1 ,r 2 (η r 1 , η r 2 ) = 1 ν ν ℓ=1 exp{- 1 2 s 2 n 2 k=1 ( s n s n η exp,ℓ r k -η r k ) 2 } . (3.41)
In (3.39) and (3.41), the constants s n and s n are given by (2.14).

Reformulation of the optimization problem.

Let i ref ≥ 0 be any fixed real value of the mutual information for two real-valued random variables. Let r 1 and r 2 be any two integers such that 1 ≤ r 1 < r 2 ≤ n. The random variables H ν r 1 and H ν r 2 will be considered as

dependent if 0 ≤ i ref < i ν (H ν r 1 , H ν r 2 ) , (3.42) independent if 0 ≤ i ν (H ν r 1 , H ν r 2 ) ≤ i ref .
We then introduce the following rule given by the probability theory: If r 1 , r 2 , and r 3 are three indices such that

1 ≤ r 1 < r 2 < r 3 ≤ n, then • if H ν r 1 is dependent on H ν r 2 and if H ν r 1 is dependent on H ν r 3 , then H ν r 1 , H ν r 2
, and H ν r 3 are mutually dependent;

• if H ν r 1 is dependent on H ν r 2 and if H ν r 2 is dependent on H ν r 3 , then H ν r 1 , H ν r 2
, and H ν r 3 are mutually dependent.

From (3.42) and this rule, it can be deduced that,

• when

i ref = 0, if 0 < i ν (H ν r 1 , H ν r 2 )
for all 1 ≤ r 1 < r 2 < r 3 ≤ n, then all the components H ν 1 , . . . , H ν n of H ν are considered as mutually dependent and, consequently, m = 1 (there is only one subset in the partition, which coincides with H ν );

• when i ref > max 1≤r 1 <r 2 ≤n i ν (H ν r 1 , H ν r 2 ), if i ν (H ν r 1 , H ν r 2 ) ≤ i ref for all 1 ≤ r 1 < r 2 < r 3 ≤
n, then all the components H ν 1 , . . . , H ν n of H ν are considered as independent and consequently, m = n (there are n independent subsets in the partition);

• when i ref is such that 0 < i ref < max 1≤r 1 <r 2 ≤n i ν (H ν r 1 , H ν r 2 )
, the construction of the partition P ν i ref (m; µ 1 , . . . , µ m ) of H ν in terms of random vectors Y ν,1 , . . . , Y ν,m considered as independent for the given level i ref , is deduced from the previous rule and from graph theory [START_REF] Bondy | Graph Theory with Applications[END_REF]. The set of nodes of the graph is I = {1, . . . , n}. Let [N ] be the symmetric adjacency (n × n) matrix made up of 0 and 1, for which all its diagonal entries are 0 and such that, for all integers r 1 and r 2 verifying 1 ≤ r 1 < r 2 ≤ n,

[N ] r 1 r 2 = [N ] r 2 r 1 = 0 if i ν (H ν r 1 , H ν r 2 ) ≤ i ref , (3.43) [N ] r 1 r 2 = [N ] r 2 r 1 = 1 if i ref < i ν (H ν r 1 , H ν r 2 )
. For such a given level i ref and from the knowledge of adjacency matrix [N ], the algorithm presented in Section 3.6 allows for constructing the unique partition P ν i ref (m; µ 1 , . . . , µ m ) associated with i ref and then for computing the numerical criterion τ ν (m; µ 1 , . . . , µ m ), defined by (3.32), which is rewritten as τ ν (i ref ),

τ ν (i ref ) = 1 - i ν (Y ν,1 , . . . , Y ν,m ) i ν (G 1 , . . . , G m ) . ( 3
i opt ref = inf i ref { arg max i ref ≥0 τ ν (i ref )} . (3.46)
Remarks.

• In (3.46), the inf operator is introduced for the same reason as the one given for justifying (3.33). Let T = arg max τ ν (i ref ) be the set of all the i ref ≥ 0 such that τ ν (i ref ) is maximum.

The numerical experiments have shown that set T is not always reduced to a single point. The solution is thus chosen as

i opt ref = inf{T }. • If for all i ref > 0, we have τ ν (i ref ) ≤ τ (i opt ref ) < 0 with |τ (i opt ref )| ≤ ǫ ≪ 1 (for instance, ǫ = 10 -3
), then all the components H ν 1 , . . . , H ν n of H ν are Gaussian and independent. • The algorithm used for solving the optimization problem defined by (3.46) is the trial method.

• Equation (3.45) means that, once the optimal value i opt ref of the level i ref has been estimated by (3.46), then the partition

P ν i opt ref (m; µ 1 , . . . , µ m ) of H ν in terms of random vectors Y ν,1 , . . . , Y ν,m
that are assumed to be mutually independent for this optimal value i opt ref of i ref , is calculated and, consequently, the optimal values m max and µ opt 1 , . . . , µ opt mmax such that µ opt 1 + . . . + µ opt mmax = n are deduced.

Algorithm for constructing the partition

P ν i ref (m; µ 1 , . . . , µ m
) of H ν in terms of random vectors considered as independent for a given level i ref .

For i ref such that 0 < i ref < max 1≤r 1 <r 2 ≤n i ν (H ν r 1 , H ν r 2
), an algorithm using the graph theory is proposed for constructing the partition P ν i ref (m; µ 1 , . . . , µ m ) of H ν in terms of random vectors Y ν,1 , . . . , Y ν,m considered as independent for the given level i ref . This algorithm, written using the Matlab language, is detailed in Algorithm 1 given in Appendix A, for which it is assumed that the symmetric adjacency (n × n) matrix [N ] has previously been calculated using (3.43). The proposed algorithm uses notions of graph theory, such as the "nodes" and the "symmetric adjacency" for which the definition can be found in [START_REF] Bondy | Graph Theory with Applications[END_REF]. The algorithm that is detailed in Appendix A has especially been developed for constructing the partition by following the development presented in Section 3.5.3. The CPU time generated by the proposed algorithm is totally due to the construction of the symmetric adjacency matrix [N ], which requires the calculation of (3.35) to (3.41) (use of the classical Gaussian kernel density estimation method for computing a realization of the random mutual information). Consequently, the CPU time is proportional to n × ν × n(n + 1)/2 (dimension of the random vector) × (number of realizations) × (number of entries in the symmetric adjacency matrix) and to the number of constructions of matrix [N ]. The CPU time generated by the graph algorithm detailed in Appendix A is completely negligible. It should be noted that the sequences that generate the CPU time can easily be parallelized due to the algebraic structure of (3.39) and (3.41), and due to the fact that the computation of the entries of matrix [N ] can also easily be parallelized.

Numerical experiments and numerical validation.

Four numerical experiments are presented. For the first three, the construction of each numerical experiment is performed in order that the optimal partition in terms of mutually independent random vectors is known a priori. Therefore, these numerical experiments yield a validation of the method proposed. For the last numerical experiment (the fourth), the optimal partition is unknown. Consequently, this numerical experiment does not directly contribute to the validation, but this example that is relative to a non-Gaussian M + 6 -valued random field indexed by R 3 , allows for testing the capability of the method proposed to identify an optimal partition for a case in high dimension. (i) Generation of the set of realizations. This numerical experiment consists in analyzing a very simple case of "almost" independent uniform real-valued random variables.

• For N = 10, let W = (W 1 , . . . , W N ) be the random vector, defined on the probability space (Θ, T , P), whose components are written as W = √ 3 (2 U -1) in which the components U 1 , . . . , U N of random vector U are N mutually independent uniform random variables on [0 , 1]. For ν = 100, let {U(θ ℓ ), ℓ = 1, . . . , ν} be ν independent realizations of U with

θ ℓ ∈ Θ. Let W(θ ℓ ) = √ 3 (2 U(θ ℓ ) -1)
be the corresponding realizations of W. • We now defined the R N -valued random vector X ν from the set of its ν independent realizations {x exp,ℓ , ℓ = 1, . . . , ν} as explained in Section 2.1, such that x exp,ℓ = W(θ ℓ ). • The reduced-order statistical model defined in Section 2.2 and the probabilistic construction of the random vector H ν given in Section 2.3 are applied in taking n = N (thus err(n, ν) = 0), which allows for constructing the ν independent realizations {η exp,ℓ , ℓ = 1, . . . , ν} of the R n -valued random vector H ν such that (2.15) and (2.16) are verified. It should be noted that, if ν went to +∞, then the components of H ν would be mutually independent. However, since the value 100 for ν is not sufficiently large, the components of H ν are not exactly mutually independent. The test for the validation of the method that is proposed is to evaluate its capability to find the mutual independence of all the components of H ν for ν = 100. As an illustration, Figure 4.1 (left) shows the graph of the pdf of the random variable H ν 3 , computed by using (3.38) with (3.39), which is compared to a Gaussian pdf with the same empirical mean value and the same standard deviation.

= (r 1 , r 2 ) → i ν (H ν r 1 , H ν r 2 )
defined by (3.35) and the horizontal line corresponding to the ordinate i opt ref = 0.02. The notation p = (r 1 , r 2 ) means that p is the integer belonging to {1, . . . , n(n -1)/2}, which is associated with the pair (r 1 , r 2 ) when r 1 and r 2 run through the ensemble {1 ≤ r 1 < r 2 ≤ n}.

The multi-indices (r 1 , r 2 ) are ordered such that for r 1 given by increasing value in 1 ≤ r 1 ≤ n -1, r 2 is ordered by increasing value such that r 1 + 1 ≤ r 2 ≤ n. This order is also used in the similar figures later in the paper. For instance, in Figure 4.2 (right), p = 1 corresponds to r 1 = 1 and r 2 = 2; p = 9, to r 1 = 1 and r 2 = 10; p = 10 to r 1 = 2 and r 2 = 3, etc. It can be seen that all the values of i ν (H ν r 1 , H ν r 2 ) are smaller than i opt ref and therefore, taking into account (3.42), there are no dependent random variables H ν r 1 with H ν r 2 for all 1 ≤ r 1 < r 2 ≤ n. We then obtained an optimal partition P ν opt = P(m max ; µ opt 1 , . . . , µ opt mmax ) with m max = n = N = 10 and µ opt 1 = . . . = µ opt mmax = 1. Consequently, the optimal partition identified corresponds to n = 10 mutually independent real-valued random variables Y ν,1 , . . . , Y ν,m with Y ν,j = H ν j , that is to say, all the components of random vector H ν are found as mutually independent for ν = 100 (that was expected). For this optimal partition, the criterion defined by (3.30) is verified with i ν (Y ν,1 , . . . , Y ν,mmax ) = 4.83 and i ν (G 1 , . . . , G m max ) = 5.29. The total CPU time for obtaining the results with a laptop computer is 7 seconds.

Numerical experiment 2.

For this numerical experiment, we test the capability of the method proposed to identify an optimal partition of a random vector having 78 non-Gaussian components, in terms of 12 mutually independent non-Gaussian random vectors of different lengths.

(i) Generation of the set of realizations. For this numerical example, we consider a random vector with N = 78 components, which is constituted of m = 12 independent random vectors with length µ j = j for j = 1, . . . , m. The generation is performed as follows.

• Let Ξ 1 , . . . , Ξ m be m mutually independent real-valued random variables (the random germs) such that Ξ j = (2 U j -1)/4 in which U 1 , . . . , U m are m mutually independent uniform random variables on [0 , 1], defined on the probability space (Θ, T , P). The number of independent realizations in the set that we are constructing is ν = 100. For ℓ = 1, . . . , ν, the ℓ-th independent realization of {Ξ 1 , . . . , Ξ m } is denoted by {ξ ℓ 1 , . . . , ξ ℓ m }. • For j fixed in {1, . . . , m}, we define the random vector Z j = (Z j 1 , . . . , Z j µ j ) with µ j = j such that, for all k = 1, . . . , µ j ,

Z j k = k i=0 √ i! (Ξ j ) i . Let m j Z ∈ R µ j and [ C j Z ]
∈ M µ j be the empirical estimates of the mean value and of the covariance matrix of Z j , computed by using the ν independent realizations ξ 1 j , . . . , ξ ν j of Ξ j . Let [L j Z ] be the upper triangular matrix such that [

C j Z ] = [L j Z ] T [L j Z ]. • We introduce the normalized R µ j -random vector W j such that W j = [L j Z ] -T (Z j -m j Z
) for which we compute the ν independent realizations {w j,ℓ ∈ R µ j , ℓ = 1, . . . , ν}. The empirical estimates of the mean value m j W ∈ R µ j and of the covariance matrix [ C j W ] ∈ M µ j of random vector W j , computed by using {w j,ℓ , ℓ = 1, . . . , ν} are such that m j W = 0 µ j and [

C j W ] = [I µ j ].
• We now introduce the R N -random vector W = (W 1 , . . . , W m ) and we denote by {w ℓ , ℓ = 1, . . . , ν} its ν independent realizations w ℓ = (w 1,ℓ , . . . , w m,ℓ ) ∈ R N . • The R N -valued random vector X ν is now defined by its ν independent realizations {x exp,ℓ , ℓ = 1, . . . , ν} as explained in Section 2.1 such that x exp,ℓ = w ℓ . • Let m ν X ∈ R N and [ C ν X ] ∈ M N be the empirical estimates of the mean value and of the covariance matrix of X ν , computed by using {x exp,ℓ , ℓ = 1, . . . , ν} (see (2.1)). We thus have m ν X = 0 N , but although [ C j W ] = [I µ j ] for all j = 1, . . . , m, since the value 100 of ν is not very big, the statistical estimate [ C ν X ] of the covariance matrix of X ν is not equal to the identity matrix and we have only

[ C ν X ] ≃ [I N ] (the diagonal blocks are effectively [I µ 1 ], . . . , [I µm ],
but the entries of the extra-diagonal blocks are not exactly zero).

• Despite the fact that [ C ν X ] is not exactly equal to [I N ], in order to preserve the mutual independence for the m random vectors that we are generating, the reduced-order statistical model defined in Section 2.2 is not applied and we define the random vector H ν such that H ν = X ν with n = N , and therefore, the ν independent realizations {η exp,ℓ ∈ R n , ℓ = 1, . . . , ν} of H ν are such that such that (2.15) is verified but (2.16) is approximatively verified.

As an illustration, Figure 4.3 shows the graph of the pdf for two components of H ν , computed by using (3.38) with (3.39), which is compared to a Gaussian pdf with the same empirical mean value and with the same standard deviation. Figure 4.4 displays the graph of the joint pdf of the random variables (H ν r 1 , H ν r 2 ) for two pairs of indices, computed by using (3.40) with (3.41). Note that the first one (left figure) corresponds to two dependent components belonging to subset j = 6 while the second pair corresponds to two independent components that belong to subsets j = 9 and j = 11 (see (4.1)). Therefore, the construction proposed allows for writing the following partition H ν = (Y ν,1 , . . . , Y ν,m ) in which Y ν,j is a random vector with values in R µ j with µ j = j. The ν independent realizations of Y ν,j are {η j,exp,ℓ , ℓ = 1, . . . , ν} such that η exp,ℓ = (η 1,exp,ℓ , . . . , η m,exp,ℓ ) with η j,exp,ℓ = (η j,exp,ℓ 1 , . . . , η j,exp,ℓ µ j ) = (η exp,ℓ r j 1 , . . . , η exp,ℓ r j µ j

) ∈ R µ j . With the construction proposed, if ν was sufficiently big, then the number of mutually independent random vectors would be surely m = 12 and, for j = 1, . . . , m, the number of components of Y ν,j = (H ν r j 1 , . . . , H ν r j µ j

) would be surely µ j = j. Let {[S] jk } jk be the (m × m) matrix of integers for which only the entries 1 ≤ j ≤ k ≤ m are defined and such that [S] jk = r j k where r j 1 , . . . , r j µ j are the indices of the mutually dependent components that are the components of Y ν,j ,

[S] =                      1 2
                     . ( 4.1)
For instance, for j = 4, the indices of H ν that constitute the components of Y ). As previously commented, It should be noted that, if ν was sufficiently big, then we should find as the optimal partition of H ν , the one defined by (4.1). Nevertheless, as the value 100 for ν is not sufficiently big, the random vectors Y ν,1 , . . . , Y ν,m defined above are not exactly mutually independent. In the framework of the validation of the method proposed, this numerical experiment consists in evaluating the capability of the method to identify the optimal partition defined by (4.1) with a small value of ν (ν = 100). Note that the [S] matrix obtained with the optimal partitioning using the proposed algorithm is exactly the one defined by (4.1).

= (r 1 , r 2 ) → i ν (H ν r 1 , H ν r 2 )
defined by (3.35) and the horizontal line corresponds to the ordinate i opt ref = 0.07. The optimal partition P ν opt = P(m max ; µ opt 1 , . . . , µ opt mmax ) that is identified with the method proposed yields m max = 12 and µ opt j = j with j = 1, . . . , 12. Consequently, the partition defined by (4.1) is identified by the method proposed for ν = 100 (as expected). For this optimal partition, the criterion defined by (3.30) is verified with i ν (Y ν,1 , . . . , Y ν,mmax ) = 23.87 and i ν (G 1 , . . . , G m max ) = 37.02. The total CPU time for obtaining the results with a laptop computer is 105 seconds. 4.3. Numerical experiment 3. This numerical experiment is the same as numerical experiment 2, except that the uniform random germ is replaced by a Gaussian random germ, but, obviously, random vector H ν is not Gaussian. This experiment is carried out in order to analyze the sensitivity of the method proposed with respect to the probability distribution of the random vector defined by the set of realizations. The reader will see that the graph of the function

i ref → τ ν (i ref ), the graph of the function p = (r 1 , r 2 ) → i ν (H ν r 1 , H ν r 2 )
, and the value i opt ref are different from the results obtained with the numerical experiment 2 (that is normal because the probability distribution is not the same) but, however, the method proposed gives the solution, which is the same as the one for numerical experiment 2. (i) Generation of the set of realizations. For this numerical example, we consider the generation presented in Section 4.2 for which Ξ 1 , . . . , Ξ m are m mutually independent real-valued normalized Gaussian random variables (centered and unit variance). All the other parameters are those of numerical experiment 2. As an illustration, Figure 4.6 displays the graph of the pdf for two components of H ν , computed by using (3.38) with (3.39), which is compared to a Gaussian pdf with the same empirical mean value and with the same standard deviation. p = (r 1 , r 2 ) → i ν (H ν r 1 , H ν r 2 ) defined by (3.35) and the horizontal line corresponding to the ordinate i opt ref = 0.13. The optimal partition P ν opt = P(m max ; µ opt 1 , . . . , µ opt mmax ) that is identified with the method proposed yields m max = 12 and µ opt j = j with j = 1, . . . , 12. Consequently, the partition defined by (4.1) is identified by the method proposed for ν = 100 (as expected). For this optimal partition, the criterion defined by (3.30) is verified with i ν (Y ν,1 , . . . , b f Y ν,mmax ) = 19.55 and i ν (G 1 , . . . , G m max ) = 37.02. The total CPU time for obtaining the results with a laptop computer is 140 seconds.

Numerical experiment 4.

For this numerical experiment, there are no references that allow for validating the optimal partition identified by the method proposed, but the example is in high dimension. The data set, which defines the random vector X ν that has N = 27, 951 components, is constituted of ν = 1, 000 independent realizations {x exp,ℓ , ℓ = 1, . . . , ν} with x exp,ℓ ∈ R N . This random vector is generated as the spatial discretization of a non-Gaussian positive-definite matrix-valued random field. ), the number m max of mutually independent random vectors for the optimal partition P ν opt = P(m max ; µ opt 1 , . . . , µ opt mmax ), the rate r of mutually independent random vectors defined as r = m max /n, and for the optimal partition, the values of i ν (Y ν,1 , . . . , Y ν,mmax ) and i ν (G 1 , . . . , G mmax ) that are used for the criterion defined by (3.30). It can be seen that rate r of independent random vectors increases when the correlation length L corr decreases, which is physically coherent; • the optimal partition defined by Figure 4.12 that displays the graph j → {[S] jk = r j k , 1 ≤ k ≤ µ opt j } in which r j 1 , . . . , r j For L corr = 0.2, for which n = 706 and ν = 1000, the CPU time with one core for constructing one symmetric adjacency matrix is about 500 minutes, but as explained in Section 3.6, with J cores, the elapsed time would be 500/J minutes (for the presented applications, the computation has be carried out with J = 12).

The results presented in table 4.1 and in figure 4.12 show that a very important gain would be obtained for the construction of a stochastic representation of such a non-Gaussian tensor-valued random field for which a set of its realizations would be given. For instance, instead of directly identifying the coefficients of its PCE from the set {η exp,ℓ ∈ R n , ℓ = 1, . . . , ν} of its realizations, this preprocessing of the data set allows for identifying, independently, for each j fixed in {1, . . . , m}, the coefficients of the polynomial chaos expansion of the random vector Y ν,j with values in R µ j for which µ j < n, using solely its realizations {η j,exp,ℓ ∈ R µ j , ℓ = 1, . . . , ν} that are directly extracted from the set {η exp,ℓ ∈ R n , ℓ = 1, . . . , ν} as soon as the optimal partition of H ν has been identified with the algorithm presented. Consequently, the methods proposed, for instance in [START_REF] Arnst | Identification of Bayesian posteriors for coefficients of chaos expansion[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data[END_REF], could be reused with a larger efficiency. Note that the partition that is constructed is relative to vector H ν , for which the components are non-Gaussian, centered, and uncorrelated (but statistically dependent) random variables that are related to the representation of the M ) and which is a strongly non-Gaussian M + 6 -valued random field (see [START_REF] Soize | Non Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF]). This means that, even if all the uncorrelated random components of H ν were mutually independent (and consequently, H ν would be a Gaussian random vector), the components of the non-Gaussian random vector X ν related to the M + 6 -valued random field [G] would be dependent. As shown in Figure 4.9 (central and right) and in Figure 4.10, some components of random vector H ν are relatively close to Gaussian random variables but are not Gaussian. This is a reasonable explanation for justifying the relatively large values of parameter r given in Table 4.1, but again, it should be noted the coherence of the results obtained, because r decreases when the spatial correlation length of the random field increases.

Conclusions.

In this paper, we have proposed a fast algorithm based on the use of mutual entropy from information theory and on the use of graph theory for constructing an optimal partition, in terms of mutually independent random vectors, of the components of a non-Gaussian random vector that is only defined by a given set of realizations. The method has especially been developed for random vectors in high dimension and for which the number of realizations that constitute the data set can be small. The method proposed and its objective are different from the ICA that was introduced to extract independent source signals from a linear mixture of independent stochastic processes, which is used for source-separation problems. The statistical tool that is proposed allows for improving the identification of any stochastic model of a random vector in high dimension for which a data set of realizations is given. Instead of directly constructing a unique stochastic model for which its stochastic dimension, which is identified by solving a statistical inverse problem, can be large, the proposed preprocessing of the data set allows for constructing several mutually independent stochastic models with smaller stochastic dimensions. Consequently, such a method allows for decreasing the cost of the identification and/or to make possible an identification for a case that is a priori in high dimension and that could not be identified through a direct and global approach. For instance, instead of directly identifying the coefficients of a PCE for a random vector from a given set of its realizations, the preprocessing of the data set allows for identifying the coefficients of several PCEs for several random vectors that have smaller dimensions. The random vector defined by the given data set can come from the discretization of a stochastic process or of a random field for which the set of realizations come from experimental measurements or computational simulations. The algorithm is completely defined in the paper, which can easily be implemented. Three numerical experiments have been proposed to validate the proposed approach. Clearly, many other cases could be analyzed. A non-Gaussian matrix-valued random field has been analyzed and the results obtained show that a very important gain can be obtained for the construction of a stochastic representation of this random field. Finally, it should be noted that the sole information that is supposed to be available for the random vector considered, is a limited set of its realizations that come, for example, from experimental measurements, and for which one seeks to build a stochastic representation by solving for example, a statistical inverse problem. It is assumed that no further information is available. So if a stochastic model is identified to represent this random vector, it will depend solely on this limited number of realizations. It is therefore consistent, in such a construction process, to use the same limited set of data to separate its components as an ensemble of several mutually independent random vectors that have smaller dimensions, before discussing the construction of a stochastic model for its representation.

  [x]. [x] T : transpose of matrix [x]. tr{[x]}: trace of a square matrix [x]. [I n ]: identity matrix in M n .

  ) in which the conventions 0 log(0/a) = 0 for a ≥ 0 and b log(b/0) = +∞ for b > 0 are used. The mutual information is related to the relative entropy (Kullback-Leibler divergence) of the pdf p Y ν,1 ,...,Y ν,m with respect to the pdf p

(

  ii) From (3.11) and (3.12), it can be deduced that S(p Y ν,1 ,...,Y ν,m ) ≤ S(p Y ν,1 ) + . . . + S(p Y ν,m ) . (3.13) (iii) If random vectors Y ν,1 , . . . , Y ν,m are mutually independent, then from (3.4) and (3.9), it can be deduced that i(Y ν,1 , . . . , Y ν,m ) = 0 . (3.14) Therefore, (3.11) and (3.14) yield S(p Y ν,1 ,...,Y ν,m ) = S(p Y ν,1 ) + . . . + S(p Y ν,m ) . (3.15) 3.3.4. Defining a theoretical criterion for testing the mutual independence of the random vectors of a partition. Let us assume that n and ν are fixed such that (2.7) is satisfied.

. 27 )

 27 Equation(3.16) shows that E{Z ν } = 1 and (3.24) shows that B ν tends tolog(1). Consequently, lim ν→+∞ B ν = 0 in probability .(3.28)

  ν in terms of mutually independent random vectors. Using the numerical criterion defined by (3.32), we propose to construct the solution of the problem formulated in Section 3.2 by solving the following optimization problem: (m max ; µ opt 1 , . . . , µ opt mmax ) = inf m {arg max 1≤m≤n { max (µ 1 ,...,µm)∈Cm τ ν (m; µ 1 , . . . , µ m )}} , (3.33) in which the admissible set C m ⊂ {1, . . . , n} m is such that C m = {µ 1 ≥ 1, . . . , µ m ≥ 1 , µ 1 + . . . + µ m = n} and where the notation used in (3.33) is detailed in the fourth paragraph of the following remarks.

4. 1 .

 1 Numerical experiment 1. The generation of the set of realizations and the results are presented below.

  Figure 4.1 (right) displays the graph of the joint pdf of the random variables (H ν 1 , H ν 3 ) computed by using (3.40) with (3.41). (ii) Results. The graph of function i ref → τ ν (i ref ), defined by (3.44) is displayed in Figure 4.2 (left). The optimal level i opt ref of i ref , defined by (3.46), is i opt ref = 0.02.

Figure 4 . 2 (

 42 right) displays the graph p

Figure 4 . 1 .

 41 Figure 4.1. Numerical experiment 1. Graph of η3 → p H ν 3 (η3) (thick line) compared to the graph of a Gaussian pdf with the same empirical mean value and with the same standard deviation as H ν 3 (thin line) (left figure) and graph of (η1, η3) → p H ν 1 ,H ν 3 (η1, η3) (right figure).

1 ,H ν r 2 )Figure 4 . 2 .

 1242 Figure 4.2. Numerical experiment 1. Graph of function i ref → τν (i ref ) (left figure) and graph of p = (r1, r2) → i ν (H ν r 1 , H ν r 2 ) (symbols) and value of i opt ref (horizontal line) (right figure).

Figure 4 . 3 .

 43 Figure 4.3. Numerical experiment 2. Graph of ηj → p H ν j (ηj ) (thick line) compared to the graph of a Gaussian pdf with the same empirical mean value and with the same standard deviation as H ν j (thin line), for j = 17 (left figure) and j = 57 (right figure).

Figure 4 . 4 . 1 ,

 441 Figure 4.4. Numerical experiment 2. Graph of the joint pdf (ηr 1 , ηr 2 ) → pHν r 1 ,H ν r 2 (ηr 1 , ηr 2 ) for r1 = 17 and r2 = 19 (left figure), r1 = 38 and r2 = 56 (right figure).

Figure 4 . 5 .

 45 Figure 4.5. Numerical experiment 2. Graph of function i ref → τν (i ref ) (left figure) and graph of p = (r1, r2) → i ν (H ν r 1 , H ν r 2 ) (symbols) and value of i opt ref (horizontal line) (right figure).

Figure 4 . 5 (

 45 right) displays the graph p

Figure 4 . 6 .

 46 Figure 4.6. Numerical experiment 3. Graph of ηj → p H ν j (ηj ) (thick line) compared to the graph of a Gaussian pdf with the same empirical mean value and with the same standard deviation as H ν j (thin line), for j = 57 (left figure) and j = 59 (right figure).

Figure 4 .

 4 7 displays the graph of the joint pdf of the random variables (H ν r 1 , H ν r 2 ) for two pairs of indices, computed by using (3.40) with(3.41). Note that the first one (left figure) corresponds to two dependent components belonging to subset j = 11 while the second pair corresponds to two independent components that belong to subsets j = 5 and j = 8 (see (4.1)).

(

  ii) Results. The graph of function i ref → τ ν (i ref ), defined by (3.44), is displayed in Figure 4.8 (left). The optimal level i opt ref of i ref , defined by (3.46), is i opt ref = 0.13.

Figure 4 Figure 4 . 7 .

 447 Figure 4.7. Numerical experiment 3. Graph of the joint pdf (ηr 1 , ηr 2 ) → pHν r 1 ,H ν r 2 (ηr 1 , ηr 2 ) for r1 = 57 and r2 = 59 (left figure), r1 = 11 and r2 = 29 (right figure).

1 ,H ν r 2 )Figure 4 . 8 .

 1248 Figure 4.8. Numerical experiment 3. Graph of function i ref → τν (i ref ) (left figure) and graph of p = (r1, r2) → i ν (H ν r 1 , H ν r 2 ) (symbols) and value of i opt ref (horizontal line) (right figure).

Figure 4 . 9 .

 49 Figure 4.9. Numerical experiment 4. Graph of the relative error function n → err(n, ν) defined by (2.6) for the three values of Lcorr (left figure): 0.2 (blue thin line), 0.5 (black thin line), and 1.0 (thick blue line). For Lcorr = 1.0, graph of ηj → p H ν j (ηj ) (thick line) compared to the graph of a Gaussian pdf with the same empirical mean value and with the same standard deviation as H ν j (thin line), for j = 33 (central figure) and j = 48 (right figure).

Figure 4 . 10 .

 410 Figure 4.10. Numerical experiment 4. For Lcorr = 1.0, graph of the joint pdf (ηr 1 , ηr 2 ) → pHν r 1 ,H ν r 2 (ηr 1 , ηr 2 ) for r1 = 17 and r2 = 33 (left figure) and r1 = 48 and r2 = 56 (right figure).

µ opt j are the indices of the mutually dependent components H ν r j 1 ,

 1 H ν , which belong to the subset j (see(3.2)).

Figure 4 . 11 .

 411 Figure 4.11. Numerical experiment 4. Graph of function i ref → τν (i ref ) (left figure) for Lcorr = 0.2 (blue thin line), = 0.5 (black thin line), = 1.0 (blue thick line) and, for Lcorr = 1.0, graph of p = (r1, r2) → i ν (H ν r 1 , H ν r 2 ) (symbols) and value of i opt ref (horizontal line) (right figure).

S 6 -Figure 4 . 12 .

 6412 Figure 4.12. Numerical experiment 4. Graph of j → {[S] jk = r j k , 1 ≤ k ≤ µ opt j } corresponding to the optimal partition for Lcorr = 0.2 (left figure), Lcorr = 0.5 (central figure), Lcorr = 1.0 (right figure).

  .44) Consequently, the optimization problem defined by (3.33) can be rewritten as the following one, P ν (m max ; µ opt 1 , . . . , µ opt mmax ) = P ν

	i	ref opt	(m; µ 1 , . . . , µ m ) ,	(3.45)
	in which the optimal level i opt ref of i ref is given by

  4,ν are 7, 8, 9, and 10, which means that Y 4,ν = (H ν 7 , H ν 8 , H ν 9 , H ν 10

Table 4 . 1

 41 Results obtained for numerical experiment 4

	L corr	n	i opt ref	m max	r	i ν (Y ν,1 , . . . , Y ν,mmax ) i ν (G 1 , . . . , G mmax )
	0.2 706 0.0084 563 0.7975	521.04	521.71
	0.5 245 0.0082 187 0.7633	183.54	183.86
	1.0	77 0.0095	51	0.6623	56.05	56.31
	•

Table 4 .

 4 1 yields the dimension n of H ν , the optimal level i opt ref of i ref , defined by(3.46

(i) Generation of the set of realizations. We begin by defining the stochastic model of the random field and then we explain how the set of realizations is constructed.

• Let {[G(ζ)], ζ ∈ R 3 } be the non-Gaussian second-order random field, defined on the probability space (Θ, T , P), indexed by R 3 , with values in M + 6 , which is completely defined in [START_REF] Soize | Non Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF] (in this reference, this field [G] is denoted as [G n ]). The stochastic model of this field depends on hyperparameters that are (1) a dispersion parameter δ that allows for controlling the level of statistical fluctuations of the random field, and (2) three lengths L 1 , L 2 , and L 3 that allow for controlling the three spatial correlation lengths of the random field, associated with the three coordinates ζ 1 , ζ 2 , ζ 3 of ζ. In the present numerical experiment, δ = 0.7 and L 1 = L 2 = L 3 = L corr for which three values L corr = 0.2, 0.5, and 1.0 will be considered.

The random generator presented in [START_REF] Soize | Non Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF] allows for constructing the ν = 1, 000 independent realizations {g ℓ , ℓ = 1, . . . , ν} with g ℓ ∈ R N of the non-Gaussian random vector G log resh .

• Let m ν G be the empirical mean value of G log resh computed by using {g ℓ , ℓ = 1, . . . , ν}. We then introduced the ν centered realizations {x exp,ℓ , ℓ = 1, . . . , ν} such that x exp,ℓ = g ℓm ν G . We now defined the R N -valued random vector X ν from the set of the ν independent realizations {x exp,ℓ , ℓ = 1, . . . , ν} as explained in Section 2.1, such that its independent realizations are X ν (θ ℓ ) = x exp,ℓ . Note that the empirical mean value of X ν is 0 N . • For each value of L corr , a reduced-order statistical model X (n,ν) of X ν is carried out as explained in Section 2.2. Taking into account the dimension N , the empirical covariance matrix [ C ν X ] defined by (2.1) with m ν X = 0 N is not computed and the reduced-statistical model is constructed by using the method derived from [START_REF] Andrews | Singular value decomposition and digital image processing[END_REF][START_REF] Gerbrands | On the relationships between SVD, KLT and PCA[END_REF][START_REF] Huang | Picture Processing and Digital Filtering[END_REF], which is based on a "thin SVD" [START_REF] Golub | Matrix Computations[END_REF] of the rectangular (N × ν) matrix [x exp,1 . . . x exp,ν ]. For the three values of L corr , Figure 4.9 (left) displays the graph of the function n → err(n, ν) defined by (2.6). The value n of the reduced-order statistical model is chosen in order that the relative error err(n, ν) is equal to 0.99, which yields n = 706, 245, and 77 for L corr = 0.2, 0.5, and 1.0, respectively. As an illustration, for L corr = 1.0, Figure 4.9 (central and right) displays the graph of the pdf for two components of H ν , computed by using (3.38) with (3.39), which is compared to a Gaussian pdf with the same empirical mean value and the same standard deviation. Figure 4.10 displays the graph of the joint pdf of the random variables (H ν r 1 , H ν r 2 ) for two pairs of indices, computed by using (3.40) with (3.41). The optimal partition obtained, for which the results are presented hereinafter, shows that the first one (left figure) corresponds to two dependent components H ν 17 and H ν 33 that belong to subset j = 14 while the second pair corresponds to two independent components H ν 48 and H ν 56 that belong to subsets j = 34 and j = 3 (see figure 4.12-right).

Appendix A. Algorithm 1.

In this appendix, we give Algorithm 1 that allows for constructing the partition P ν i ref (m; µ 1 , . . . , µ m ) of H ν in terms of random vectors that are considered as independent for a given level i ref . This basic algorithm is derived from the graph theory and uses the symmetric adjacency matrix [N ]. The algorithm is described using the Matlab language.

-The input variables are: n and, for i ref given, the (n × n) symmetric adjacency matrix [N ].

-The output variables are: m, the m-column matrix M such that M (j) = µ j for j = 1, . . . , m, and the (m × max j µ j )-matrix [S] such that the j-th row of [S] is made up of the indices [r j 1 . . . r j µ j ] of the random vector

).

Algorithm 1 j = 0; % M = zeros(n, 1); S = zeros(n, n); % M (j) = µ j ; S(j, k) = r j k (see (3.2)). U = zeros(n, 1); % U (r) = 0, node r not treated and = 1, treated. while isempty(find(U = = 0)) = = 0 % do: there are nodes that have not been treated. U = find(U = = 0); % nodes that have not been treated. x = U (1); % node used for starting a new subset.

P = [ ];

% initializing the list of the nodes to be analyzed.

V = [ ];

% initializing the list of nodes already analyzed.

R = [x];

% R contains the nodes of the present subset. j = j + 1;

% increasing the present number of subsets. P = union(P, x); % union of the nodes with no repetitions. while isempty(P ) = = 0 % P not empty do. y = P (1); P (1) = [ ];

% load a node and unstack P . V = union(V, y); % union of the nodes with no repetitions. for z = 1 : n % exploring all nodes z such that N (y, z) = 1 W = find(union(P,

% y belongs to the subset. end % end if. end % end for. end % end while. µ j = size(R, 2); M (j) = µ j ; S(j, 1 : µ j ) = R; % U (R) = 1; % all nodes in R have then been treated. N (R, R) = 0; % setting to zero the nodes in subset j.