Louis Pujol
email: louis.pujol@universite-paris-saclay.fr

ISDE: Independence Structure Density Estimation

Keywords: Independence Structure Density Estimation Multivariate Density Estimation, Independence Structure, Computational Statistics

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

NOTATIONS

Let f be a density function (a nonnegative real function whose integral is equal to 1) over R d . If we think of f from a statistical viewpoint, it is natural to refer to the indices {1, . . . , d} as the features.

Let S ⊂ {1, . . . , d}, we denote by f S the marginal density of f over S. For all x = (x 1 , . . . ,

x d ) ∈ R d f S (x) = f (x) i / ∈S dx i . (1)
With a slight abuse of notation, to highlight the fact that f S (x) does not depend on (x i ) i / ∈S , we write f S (x S ) instead of f S (x).

Let k be an positive integer not greater than d. We denote by Set k d the set of all subsets of {1, . . . , d} with cardinal not greater than k and by Part k d the collection of all partitions of {1, . . . , d} constructed with blocks in Set k d . We also use the shortcuts Set d = Set

INTRODUCTION

Unsupervised Learning and Density Estimation Unsupervised learning is an important field of data analysis. It aims to design methods to extract meaningful information from a dataset with little prior knowledge. A central task in unsupervised learning is density estimation. Given a sample X 1 , . . . , X N drawn independently from a random variable X on R d with a density f , the goal is to build an estimator f of f . This question finds many applications, and density estimation is a building block for many learning tasks such as clustering ( [START_REF] Chazal | Persistence-based clustering in riemannian manifolds[END_REF], [START_REF] Campello | Density-based clustering based on hierarchical density estimates[END_REF]) or anomaly detection ( [START_REF] Chandola | Anomaly detection: A survey[END_REF]) among others.

Nonparametric and Parametric Density Estimation

The easiest way to do density estimation is to consider parametric models: data is supposed to be drawn from a probability distribution known up to a finite-dimensional parameter 1 θ. Estimating the density is then equivalent to estimating θ. One example is the centered multivariate Gaussian framework, where the parameter θ is the covariance matrix Σ. An introduction to parametric statistics can be found in [START_REF] Wasserman | All of statistics: a concise course in statistical inference[END_REF], chapter 9. This approach suffers from a lack of flexibility as it strongly constrains the model. At the other end of the spectrum lies nonparametric density estimation. In this framework, densities are no longer considered members of some finite-dimensional family but are supposed to belong to a set of functions with a given regularity. An introduction to the subject can be found in [START_REF] Alexandre | Introduction to Nonparametric Estimation[END_REF].

Kernel Density Estimators In the sequel, we focus on nonparametric density estimation. Kernel Density Estimator (KDE) is a popular density estimator in this context. It has its origins in the works of Rosenblatt [START_REF] Rosenblatt | Remarks on Some Nonparametric Estimates of a Density Function[END_REF] and Parzen [START_REF] Parzen | On estimation of a probability density function and mode[END_REF]. It has been successfully used to real-world applications in recent years (connectivity among salmon farms [START_REF] Cantrell | The use of kernel density estimation with a bio-physical model provides a method to quantify connectivity among salmon farms: spatial planning and management with epidemiological relevance[END_REF], physical activity [START_REF] King | The use of kernel density estimation to examine associations between neighborhood destination intensity and walking and physical activity[END_REF], ecological niche modelling [START_REF] Qiao | A cautionary note on the use of hypervolume kernel density estimators in ecological niche modelling[END_REF], modelling of T cell receptors [START_REF] Owens | T cell receptor signalling inspired kernel density estimation and anomaly detection[END_REF], among many others).

In this paper, we will consider Spherical Gaussian KDE (SGKDE). For a given bandwidth h > 0 we define the SGKDE associated to h and to the sample X 1 , . . . , X N as

fh (x) = 1 N N i=1 exp -(X i -x) T (X i -x) 2h 2 (2π) d/2 h d . (2)
As we will not consider other choices of kernels, we write KDE instead of SGKDE. The construction of the estimator over a data sample corresponds to the choice of the bandwidth. Different approaches exist. In practice, a cross-validation scheme over a collection of potential values of h is a popular choice. See [START_REF] Mark | Asymptotic optimality of likelihood-based cross-validation[END_REF] for analysis in the context of maximum likelihood density estimation.

Curse of Dimensionality When dealing with multidimensional data, one must be aware of the issues that the number of features can imply. It is a general fact that for the majority of statistical tasks, the higher the dimension is, the harder the estimation is (see, for example, [START_REF] Giraud | Introduction to high-dimensional statistics[END_REF]). For density estimation, the complexity can be evaluated through minimax risk, quantifying the statistical error in a worst-case scenario. It is influenced by two parameters: a regularity parameter β and dimension d , the rate of convergence for the squared L 2 loss is typically proportional to N -2β 2β+d (see [START_REF] Goldenshluger | On adaptive minimax density estimation on r d[END_REF] for a review of the literature). We remark that the higher d is, the slower the minimax risk tends to zero. This phenomenon is a manifestation of the so-called curse of dimensionality. For practitioners, it should be adventurous to use a multivariate density estimator if the sample size is limited and the dimension becomes large, especially in the case of nonparametric estimation. A solution is to assume that unknown density belongs to a class of structured functions.

Moderately High Dimension Setting

In recent years, attention was put on high-dimensional problems, where the number of features can vary from hundreds to thousands. We are interested here in situations of moderately high dimension, where the number of features can vary from a few ones to a few dozens. In this setting, the curse of dimensionality still occurs. It is of particular interest to distinguish both paradigms as we will develop algorithmic solutions that allow exhaustive search over admissible structures in moderately high dimensions but become too time-consuming in high dimensions.

Structural Density Estimation with Undirected Graphical Models A way to consider a structure for a multivariate random variable is through its undirected graphical model (introduction to the field can be found in [START_REF] Giraud | Introduction to high-dimensional statistics[END_REF] and more in-depth cover in [START_REF] Martin | Graphical models, exponential families, and variational inference[END_REF]). As we will not consider directed graphical models, we always consider that graphs are undirected in the sequel. Given a graph G = (V, E) whose vertices correspond to the features {1, . . . , d} we say that G is a graphical model for X if the following condition is satisfied:

(i, j) / ∈ E ⇒ X i ⊥ ⊥ X j |(X k ) k / ∈(i,j) . (3)
Constraints on the graphical model associated with a distribution impose a structure on the density, and such a structure can help overcome the curse of dimensionality. However, learning a graphical model is a complex task in many situations. The general result is that if G is a graphical model for a d-dimensional random variable X, denoting by C the set of cliques of G (ie fully connected sets of nodes), it exists a collection of nonnegative functions (ψ C ) C∈C such that the density f of X can be written as

f (x) = 1 Z C∈C ψ C (x c ) ( 4
)
where Z is a normalization constant. As remarked in section 2.1.2 of [START_REF] Martin | Graphical models, exponential families, and variational inference[END_REF], the functions ψ C do not have a clear relationship with the marginal densities of f . The density estimation under a graphical model for general graphs is then too ambitious, and it is necessary to constrain the graph structure.

Forest Density Estimation In a fully nonparametric setting, to our knowledge, one method is available: Forest Density Estimation (FDE) [START_REF] Liu | Forest density estimation[END_REF]. It corresponds to the estimation of a density with an uncyclic graphical model (also called a forest).

In this case, the density can be expressed with 1 and 2-dimensional marginals.

If G = (V, E) is a forest, the density f of a random variable admitting G as a graphical model enjoys the following formulation

f (x) = (i,j)∈E f {i,j} (x i , x j ) f {i} (x i )f {j} (x j ) d k=1 f {k} (x k ). (5)
In [START_REF] Liu | Forest density estimation[END_REF] the algorithm to estimate a forest and the corresponding density is presented. Let us emphasis that it requires the estimation of marginals up to dimension 2. Theorem 9 in [START_REF] Liu | Forest density estimation[END_REF] emphasis that if the true density enjoys a forest graphical model and under suitable condition on the density, the speed of convergence of FDE under Kullback-Leibler (KL) loss is related to the the speed of convergence for KDE in dimension 2 instead of in the ambient dimension d. This emphasize that FDE is a remedy to the curse of dimensionality. The KL loss between f and an estimator f is defined as

KL f f = log f f f. ( 6
)
Independence Structure In the present work, we focus on the model of Independence Structure (IS) for multivariate density developed by [START_REF] Lepski | Multivariate density estimation under sup-norm loss: oracle approach, adaptation and independence structure[END_REF] and studied by [START_REF] Rebelles | Lp adaptive estimation of an anisotropic density under independence hypothesis[END_REF]. It contains d-dimensional densities, which can be decomposed as a product of low-dimensional marginals, forming a partition of the original features.

f (x) = S∈P f S (x S ) (7)
Under a graphical model perspective, it corresponds to graphs that are composed of disjoint connected components. Previous works on IS have highlighted that if the density enjoys the property that the size of the biggest block of the partition is equal to k < d, then the complexity of density estimation, measured through minimax rate of convergence under L p losses (1 ≤ p ≤ ∞) is related to k instead of the ambient dimension d. However, these works rely on the analysis of estimators that are hardly implementable for reasonable data size.

Our Contribution We present Independence Structure Density Estimation (ISDE), a method designed to simultaneously compute a partition of the features and a density estimation as a product of marginals over this partition in order to maximize the empirical log-likelihood, or equivalently, minimize the KL loss. Our method enjoys reasonable running time for moderately high-dimensional problems and can be combined with any density estimation technique, so it covers parametric and nonparametric settings. To our knowledge, we are the first to design an algorithm estimating as IS in the context of KDE.

Organization of the Paper In section 3 we present the construction of ISDE. We compare our method with some existing ones for density estimation for synthetic datasets in section 4 and for real-world datasets in section 5 before analyzing its algorithmic complexity and running time in section 6.

ISDE

This section presents ISDE, an algorithm designed to simultaneously perform density estimation and independence partition selection in a moderately highdimensional setting.

Specifications Let k be an input parameter. We aim to provide a method taking point cloud as input and outputting an IS (a partition of the features in Part k d ) and a density estimator as a product of marginal estimators f P, ĥ P =

S∈ P fS, ĥS (8)
where ĥ P = ĥS

S k d = k i=1 d k . (9)
For B k d exact computation is harder but we can prove that (see appendix B.1) These computations indicate that it would be beneficial to find a way to avoid the computation of B k d estimators. Intuitively, as estimators are combinations of marginals estimators, it seems reasonable to decouple marginal estimations from partition selection. We will now see that we must carefully choose the loss function to implement this idea.

B k d ≥ B 2 d = 1 + d 2 + d 2 d-2 2 2! + d 2 d-2 2 d-4 2 3! • • • + d 2 . . . d-2( d/2 -

Choice of Loss Function

We have announced in the introduction that ISDE aims to minimize the Kullback-Leibler loss between the proper density and the estimate one. Here we will see that this choice is not innocuous and that other choices of loss function do not lead to a feasible algorithm.

In density estimation literature, the most popular choice for the loss function is undoubtedly the squared L 2 loss. For a partition P ∈ Part k d we want to find the collection of bandwidth ( ĥP

S ) S∈Part k d solutions of min ( ĥP S ) S∈P f -fP,h P 2 = f 2 P,h P -2 fP,h P f + f 2 . ( 11
)
If P [.] corresponds to the integral over the measure induced by the density f , an equivalent formulation is given by min

(h P S ) S∈P f 2 P,h P -2P fP,h P = min (h P S ) S∈P S∈P f 2 S,h S -2P S∈P fS,h S . (12)
Let S ∈ Set k d and P 1 , P 2 ∈ Part k d such that S ∈ P 1 and S ∈ P 2 . There is no reason to have ĥP 1 S = ĥP 2 S from the previous formulation. Then under the squared L 2 loss we have no clue on how we can avoid constructing as many estimators as elements in Part k d .

Now, for the KL loss, we want to find a collection of bandwidth ( fS ) S∈Part k d minimizing min

(h P S ) S∈P log f fP,h P f. (13)
An equivalent formulation is given by max Empirical Formulation of the Optimization Problem Under KL loss, bandwidths optimization and partition selection become two separated tasks. This decoupling incites us to design an algorithm consisting of two steps: first, compute a marginal estimator fS for all S ∈ Set k d and then find the best combination of them for a log-likelihood criterion. Let n and m be two positive integers such that m + n = N . The dataset X 1 , . . . , X N is split into two disjoint subsamples:

(
• W 1 , . . . , W m used to compute marginal estimators ( fS ) S∈Set k d • Z 1 , . . . , Z n used to compute empirical log-likelihoods ( n (S)) S∈Set k d where n (S) = 1 n n i=1 log fS (Z i )
Let us use the notation n (P) = S∈P n (S). The empirical optimization task can be written as max

P∈Part k d n (P) = max P∈Part k d S∈P n (S). ( 15
)
Partition Selection A naive approach to solve 15 is to compute n (P) for every partition of Part k d and then find the optimal one. However, this approach becomes time-consuming when d grows and infeasible for large values of d because of the number of partitions. Therefore, it will be appreciable to reformulate this optimization to speed up computation. It is possible to reformulate 15 as the following linear programming task.

Solve max

x∈R Set k d S∈Set k d n (S)x(S) (16)
Under constraints

Ax = (1, . . . , 1) T (17) x ∈ {0, 1} S k d . ( 18
)
Where x is a binary vector representing which elements of Set k d are selected, and A is a d × S k d matrix where each column is a binary vector representing the composition of one of the sets of Set k d . The condition Ax = (1, . . . , 1) T then ensures that each feature is chosen once, implying that the sets selected with x form a partition.

We validate this approach through a running time comparison (see section 3) between the implementation of a brute-force approach and a linear program solver. In this experiment, we fix the quantities ( n (S)) S∈Set k d , the brute-force approach consists in a for loop (implemented in Python), computing n (P) for all P ∈ Part k d and returning the maximum. For the LP formulation, the optimization is done with the branch-and-bound method, implemented in the Python package PulP [START_REF] Mitchell | Pulp: A linear programming toolkit for python[END_REF]. With the brute-force approach and choice k = d, partition selection takes approximately 3 hours in dimension [START_REF] Liu | Forest density estimation[END_REF]

EXPERIMENTS ON SYNTHETIC DATA

In this section, we validate the performance of ISDE on synthetic data generated under IS hypothesis.

Data Generating Process For a given list of positive integer (a structure) S = [s 1 , . . . , s K ], the data generating process is defined as follows. For each s i ∈ S, we define a s i dimensional dataset drawn from P i :

• If s i = 1, P i is the uniform distribution over [0, 1]
• If s i = 2, P i is a distribution corresponding to data sample near two concentric circles with different radii

• If s i = 3, a sample X from P i is obtained as follows: let Y 1 and Y 2 be two independent Bernoulli variables with probability of success 0.5 and

Y 3 = |Y 1 -Y 2 |. X is then drawn from the multivariate Gaussian distribution N ((Y 1 , Y 2 , Y 3 ), 0.08 × I 3 )
. This is a situation where features of P i are pairwise independent but not mutually independent

• If s i ≥ 4, P i is a mixture of two multivariate Gaussian distributions, one centered in (0, . . . , 0), the other in (1, . . . , 1)

The final dataset results from their concatenation, plus featurewise rescaling so that each value lies between 0 and 1. The dimension is d = k i=1 s i . This rescaling step does not affect the IS as it is done featurewise.

Evaluation Scheme To evaluate the performance of an estimator, we compute the empirical log-likelihood on a validation set X valid = X valid 1 , . . . , X valid M drawn independently from the same distribution as X 1 , . . . , X N :

Score( f ) = 1 M M i=1 log f X valid i . ( 19
)
The set X valid = X valid 1 , . . . , X valid M is not used to tune the estimators. In the experiments of this section, we set M = 5000.

Benchmarked Methods We will compare three density estimation algorithms for samples corresponding to different structures.

The first one is CVDKE, a KDE estimator where the bandwidth parameter is selected through a 5-fold cross-validation to maximize empirical log-likelihood on test data. The collection of possible bandwidths is a regular grid on a log-scale from 0.01 to 1 with 30 values.

The second one is ISDE with k = d (ie all partitions are tested), m = n = 0.5N and the collection of marginal estimators ( fS ) S∈Set d is a collection of CVKDE estimators constructed with the sample W 1 , . . . , W m .

The third one is FDE. Our implementation is a slight modification of the held-out data approach proposed in [START_REF] Liu | Forest density estimation[END_REF]: we rely on the quantities ( n (S)) S∈Set . We use a cross-validation scheme to optimize the bandwidth instead of the plug-in approach presented in the paper.

We insist that comparing these methods for density estimation through empirical log-likelihood for validation data is fair as all of them aim to maximize the log-likelihood.

Results Empirical log-likelihood on validation data for methods listed above are shown in table 4, for different structures and for the choice N = 5000. Each experiment is repeated 5 times, and we show the mean log-likelihood and the standard deviation on the table.

[2, 2, 1] [ For [START_REF] Campello | Density-based clustering based on hierarchical density estimates[END_REF][START_REF] Campello | Density-based clustering based on hierarchical density estimates[END_REF][START_REF] Campello | Density-based clustering based on hierarchical density estimates[END_REF], as features are pairwise independent, FDE outputs at every try a graph without any edge and computes the density as a product of one-dimensional marginals, leading to poor results in comparison to ISDE. CVKDE leads to better estimation for this setting than FDE but is outperformed by ISDE. For [4, 4, 2, 2], FDE outputs a subgraph of the actual graphical model at every try. It leads to better estimation than CVKDE but worse than ISDE, which learns the proper IS at every try. Thus, ISDE leads to better results than FDE and CVKDE for the task of structured density estimation under KL loss under IS. We interpret the bad performance of CVKDE as a manifestation of the curse of dimensionality. ISDE outperforms FDE because it considers potential higher-order dependencies between features than FDE, which only considers pairwise associations. However, let us remark that FDE covers some models not addressed by ISDE. ISDE performs better on data where IS is true, but we recommend testing both methods to determine the one that best fits the data.

We also remark that ISDE recovers exactly the IS for the considered settings. One can wonder why we do not observe that outputted partitions are not precisely the IS but partitions where blocks are a union of blocks of the true IS. We believe that this is because a useless merging of blocks in the partition is strongly penalized by ISDE as the dimension limits our ability to estimate a density accurately. Then the hold-out scheme implemented in ISDE (by splitting X into W and Z in algorithm 1) penalizes sufficiently too big blocks in partitions and leads to accurate recovery of IS.

EXPERIMENTS ON MASS CYTOMETRY DATA

This section is devoted to the presentation of some outputs on real-world datasets.

In addition to studying the performance of ISDE in terms of log-likelihood, it is the occasion to illustrate how we can interpret the outputted partition.

Datasets

The datasets presented here are the output of mass cytometry experiments. Cytometry allows high-throughput measurements at a single-cell level over a cell sample. Two types of information about cells are collected. Some are about the cell's geometry, and others about the abundance of some targeted proteins at their surface. The number of events for cytometry experiments on blood samples usually lies between 10, 000 and 1, 000, 000, and the number of features can vary from a few ones to approximately 50.

We present here results on two public cytometry datasets used in a bench-mark of clustering methods paper [START_REF] Lukas | Comparison of clustering methods for high-dimensional single-cell flow and mass cytometry data[END_REF], Levine13 and Levine32. Both are experiments on bone marrow cells extracted from healthy human donors with respectively 13 and 32 features. The preprocessing step is a featurewise rescaling to force each feature to take values in [0, 1].

Quantitative evaluation

Benchmarked Algorithms As in the previous section, we compare FDE, CVKDE, and ISDE (the value of k depends on the dimension, we selected k = 3 for Levine32 and k = 5 for Levine13 to keep computations fast).

We have also added a parametric approach to the benchmark: a Gaussian Mixture (GM) model with a selection of the number of components. This model is particularly adapted to cytometry as we naturally expect in this context that the data forms clusters representing cell populations ( [START_REF] Reiter | Clustering of cell populations in flow cytometry data using a combination of gaussian mixtures[END_REF], [START_REF] Finak | Merging mixture components for cell population identification in flow cytometry[END_REF]).

Let n C be a positive integer corresponding to the number of components in the mixture. Let p = (p 1 , . . . , p n C ) be a collection of nonnegative real number such that n C i=1 p i = 1, µ = (µ 1 , . . . , µ n C ) a collection of vector in R d and Σ = (Σ 1 , . . . , Σ n C ) a collection of d × d definite positive matrices. The density f (n C ,p,µ,Σ) of the Gaussian mixture model associated with the parameters (n C , p, µ, Σ) is

f (n C ,p,µ,Σ) = n C i=1 p i f µ i ,Σ i (20)
where f µ i ,Σ i is the density of the multivariate Gaussian random variable with mean µ i and covariance matrix Σ i .

Given n C and a dataset, it is possible to compute estimators (p, μ, Σ) with the EM algorithm [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] to maximize the log-likelihood. As we do not know the optimal number of components in advance, a strategy is to fit a Gaussian mixture model for different n C (from 1 to 30 in our experiments) and select the number of components in the mixture with a cross-validation scheme. We rely on the implementation of these methods provided by scikit-learn [START_REF] Buitinck | API design for machine learning software: experiences from the scikit-learn project[END_REF] with no restriction on the shape of the covariance matrices.

Though GM is principally used for clustering purposes, it can also be inter-preted as a parametric density estimator intended to maximize the log-likelihood. It is then relevant to compare it with the other introduced methods.

Experimental Setup From each dataset we have extracted a train sample with N = 5000 events, this train sample is exclusively used to compute estimators fCVKDE , fFDE , fISDE and fGM . For ISDE we fixed m = 3000 and n = 2000. Then to compare between these density estimators, we sampled 20 datasets with 2000 events from the data that were not used to compute estimators.

Results Boxplots indicating the log-likelihood of these estimators for validation samples can be visualized in figure 1. We remark that using ISDE leads to better empirical log-likelihood on validation data. CVKDE in the ambient dimension is always the worst estimator. GM is slightly better than FDE for both datasets, and the gap between performances of FDE/GM and ISDE is higher in dimension 32 than in dimension 13. We conclude that IS with a limited size of blocks seems to be a relevant model for these datasets as ISDE could outperform other model-based approaches in terms of log-likelihood.

Testing ISDE against other density estimation methods is a way to evaluate how this model can explain the data well. However, we must be careful in our conclusion. These results do not indicate that the data follow an IS, but rather that IS offers a good approximation of the data distribution.

Qualitative Interpretation

We believe that the added value of our method is that ISs are easy to understand and useable as a tool to interpret data. After validating the pertinence of ISDE in comparison with other methods through quantitative analysis, we now provide some insight into the capacity of ISDE to deliver meaningful qualitative information.

Nontriviality of Outputted Partition

The first question to ask is if the gain in terms of empirical log-likelihood is due to the specific outputted partition P or if any other estimator fP based on a partition of features P ∈ Part k d could achieve the same performance. To answer this question, we have computed empirical loglikelihood on 10 validation sets of size 2, 000 for the three best partitions outputted by ISDE, the three worst ones regarding the optimization task, and three random partitions in Part k d . To compute not the optimal but the second one, the third one, and so on, it suffices to add constraints on the partition selection problem that artificially exclude some partitions from the optimization. To compute the worst partitions, switching the optimization from maximization to minimization suffices. Random partitions are computed by generating a random permutation σ of {1, . . . , d} and then gather consecutive features in {σ(1), . . . , σ(d)} in groups with sizes drawn uniformly between 1 and k. These experiments indicate that ISDE outputs specific partitions that lead to better estimators in terms of log-likelihood on empirical data than the random partitions. In that sense, the information provided by ISDE on these datasets is not trivial. It also seems that not only the optimal one P but a collection of partitions lead to the best scores.

With that in mind, it could be interesting to determine if the collection of partitions leading to optimal results are close in some sense. To this end, it is necessary to introduce a notion of distance between partitions.

Edit Distance Given two partition P and P in Part k d it is possible to define a distance between P and P called edit distance ( [START_REF] Duncan P Brown | Automated protein subfamily identification and classification[END_REF]) and denoted by edit(P, P ). This distance corresponds to the minimal number of operations required to go from P to P where an operation can split a block into two or merge two blocks. The edit distance defines a distance on Part k d in the mathematical sense as it is nonnegative, symmetric, equal to zero only if we compute the distance from one partition to itself and it satisfies the triangular inequality.

Correlation between Edit Distance and Density Estimation

We will now see how the edit distance from P to P correlates with the empirical log-likelihood on validation data for fP . Secondly, we explore the space Part k d by defining a random walk considering the topology induced by edit. We define a random walk (P 0 , P 1 , . . . ) as follows: at each step we go from P i to P i+1 with edit(P i , P i+1 ) = 1. To do so, it suffices to randomly choose an operation (edit or merge) and apply it to randomly selected block(s) of P i while controlling that we stay in Part k d .

To observe a possible correlation between edit( P, .) and log-likelihood on validation data, we have implemented the following protocol: do 5 walks of length 40 with P as starting point and store all visited partitions, then for the 200 selected partitions, compute empirical log-likelihood on ten resampling of validation data and store the mean value. Then we plot these scores against edit( P, .). For both datasets, we observe a clear negative correlation between edit( P, .) and empirical log-likelihood on validation data. These observations indicate that the topology induced by the distance edit on Part k d is meaningful in the sense that the farther a partition P is from P for the edit distance, the worse the estimator fP is.

Exhaustive Analysis For the dataset Levine13, as the cardinal of Part 5 13 is 25, 719, 630, it is possible to store the entire family of empirical log-likelihood computed thanks to the data Z 1 , . . . , Z n on ISDE: ( n (P)) P∈Part 5 [START_REF] King | The use of kernel density estimation to examine associations between neighborhood destination intensity and walking and physical activity[END_REF] . Such an exhaustive analysis is impossible for Levine32 as the number of partitions in Part 3 32 exceed 10 19 . The distribution of ( n (P)) P∈Part 5

13 can be visualized thanks to an histogram. If we select the partitions with a score higher than 14.6, there remain 1, 941 elements. For these partitions, we compute empirical log-likelihood again on validation data and represent it against edit( P, .). This is a way to ask about the uniqueness of the optimal partition P. If another partition P a significantly positive value of edit( P, P) gives as good results as P, it will indicate that there are other local maximums than P. d equipped with edit distance in terms of empirical log-likelihood for fP has led us to the conclusion that the qualitative information provided by ISDE through P is nontrivial for these datasets as random partitions in Part k d does not lead to optimal scores. We also show that the density estimation score deteriorates as the edit distance from P increases, indicating that edit distance is a relevant metric to explore Part k d in density estimation under IS. Then an exhaustive analysis of the space of partitions for Levine13 indicates that we can consider the optimal partition as unique for this experiment.

These conclusions depend on the specific datasets presented here and could become invalid for other ones. We provide the code to reproduce our experiments. Our aim is that anyone interested in the method can replicate these analyses for other data.

COMPLEXITY AND RUNNING TIME ANAL-YSIS

In this section, we provide information about the algorithmic complexity and running time of ISDE.

Computation of KDE For a given bandwidth h, the evaluation of a KDE constructed over m 1 points and evaluated over m 2 points is O(m 1 m 2 ). The family of estimators fS S∈Set k d is constructed using a V -fold cross-validation where V is a divisor of m. If n h denotes the number of candidate values for the bandwidths, the number of operation required for bandwidth selection is

S k d n h V m V × m(V -1) V . The complexity of this step is O(S k d n h m 2 ).
Once the bandwidths are selected, it remains to compute the quantities ( n (S)) S∈Set k d thanks to Z 1 , . . . , Z n . The total cost of its operation is O(S k d nm). The total algorithmic cost of the computation of ( n (S))

S∈Set k d is O S k d m (n h m + n) . ( 21
)
Partition Selection The implementation of the partition selection step relies on the branch-and-bound method. It is not easy to give a precise statement about its complexity. The branch-and-bound algorithm uses a tree search strategy to enumerate all possible solutions to a given problem implicitly. A recent survey can be found in [START_REF] David R Morrison | Branch-and-bound algorithms: A survey of recent advances in searching, branching, and pruning[END_REF].

Running time We now present some information about running time. We have run all experiments on a laptop with the following hardware: CPU Intel(R) Xeon(R) W-10885M CPU @ 2.40GHz and GPU: Nvidia Quadro RTX 3000 Mobile.

The KDE computations have been performed on GPU using the python package pyKeOps [START_REF] Charlier | Kernel operations on the gpu, with autodiff, without memory overflows[END_REF]. This implementation is much faster than the one on CPU proposed by scikit learn as highlighted by table 5, which compares running time for KDE constructed on n points and evaluated on n points on dimension d = 3. The main conclusion of this running time study is that the running time of partition selection is negligible in comparison with the one for computing ( n (S)) S∈Set k d for the parameters presented here. The code associated with this paper contains functions allowing the reader to reproduce these experiments with different settings and estimate the running time on its device. Note that the code also runs if no GPU is available. In this case, pyKeOps will automatically use parallelization on CPU for KDE evaluations.

CONCLUSION

ISDE is an algorithm that outputs an estimate of a density function of a point cloud, taking into account an IS for data in moderately high dimensions. To design it, we reduced the number of hyperparameters with an appropriate choice of the loss function and, through linear programming reformulation, made the partition selection step faster than was previously possible. This leads to reasonable running time even on a laptop for the considered datasets. The code is available and ready to be used by anyone interested in this method.

ISDE is versatile: it takes any basic multidimensional density estimator as input. Then it can be used in parametric and nonparametric frameworks. It is also exhaustive as it searches over all partitions of features with given maximal block size. To our knowledge, we are the first to propose a method that considers IS in the context of nonparametric density estimation with KDE.

We validated its performance on synthetic data satisfying IS. This performance was measured in terms of log-likelihood on the validation sample. We found that ISDE exploits IS structure and outperform other density estimators for this task. Applying ISDE to mass cytometry data has indicated that it could accurately estimate density over real-world datasets and extract qualitative information about their features through the outputted partition. This paper focused on algorithmic and implementation details relative to ISDE and empirical study. Theoretical study of ISDE will be presented in a separate work, as it involves some minor modifications to prove convergence rates.

A APPENDIX: EXPERIMENTS ON GAUSSIAN SYNTHETIC DATA

This section is dedicated to the presentation of synthetic results, in the same spirit as section 4 but with data drawn from centered multivariate Gaussian distributions.

Data Generating Process

The Gaussian Graphical Models (GGM) theory indicates that edges of the undirected graphical model associated with a Gaussian distribution N (0, Σ) are the non-zero entries of the precision matrix Σ -1 . As the inverse operator preserves the block-diagonal structure, we can easily simulate data from a multivariate Gaussian with an IS.

For a positive integer s and a real number σ ∈ (0, 1) we denotes by Σ s σ the s × s matrix whose diagonal entries are 1 and nondiagonal entries are σ. Then for a list of positive integers S = [s 1 , . . . s K ] we define the block diagonal matrix:

Σ S σ =       Σ s 1 σ 0 . . . 0 0 Σ s 2 σ . . . . . . . . . . . . . . . 0 0 . . . 0 Σ s K σ       (22)
The distribution N 0, Σ S σ satisfies the IS condition with partition j-1 i=1 s i + 1, . . . , j i=1 s i j=1,...,K

.

Evaluation Scheme If Σ and Σ are respectively the estimated and the true covariance, the Kullback-Leibler risk can be explicitly computed (see appendix B.2):

KL N (0, Σ) N (0, Σ) = v∈Sp(A) v -log(1 + v) 2 (23)
where A = ( Σ-1 -Σ -1 )Σ.

Benchmarked Methods Two methods will be compared to ISDE for the task of covariance estimation.

The first estimator is the simple Empirical Covariance, which is the maximum likelihood estimator if the covariance does not enjoy any particular structure.

The second estimator is Block-Diagonal Covariance Selection (BDCS) developed in [START_REF] Devijver | Block-diagonal covariance selection for high-dimensional gaussian graphical models[END_REF]. It aims to estimate an IS in the context of GGM. This algorithm works in two steps:

• Compute a family of nested partitions candidates to be the IS • Choose a partition in this family using a slope heuristic approach More details can be found in the original paper. Up to our knowledge, this is the only work dealing specifically with IS in the GGM framework.

ISDE Inputs

We run algorithm 1 with k = d, m = n = 0.5 × N and simple empirical covariance as multivariate density estimator.

Performance We compare the three methods described above for fixed σ, N , and different structures S. We have gathered results in terms of KL loss are in table 8. We have repeated each experiment 5 times, and the scores displayed are the mean KL losses and standard deviation over these repetitions. Conclusion We remark that BDCS is the most efficient method for the task of density estimation in GGM under IS. We can explain it as ISDE tends to select admissible partition but fails to select the exact IS when the dimension grows. BDCS inherently penalizes more useless blocks merging, making it more accurate in this setting. However, ISDE performs significantly better than a naive empirical covariance, proving that it benefits from the IS.

We want to highlight the difference between ISDE and BDCS. BDCS starts by selecting a family of up to d nested partitions and then selects among them. This approach uses a preliminary covariance estimator to design this family of nested partitions. This approach is reasonable as for Gaussian data, pairwise dependencies entirely determine multidimensional dependencies between features. Outside the scope of GGM, this approach does not remain valid as features of a random variable can be pairwise independent but mutually dependent. ISDE can handle more general settings as it selects among a set of partitions with blocks of cardinal potentially more significant than 2. As selecting a partition in Part k d is equivalent to an unordered choice of pairs of variables, it remains to divide by the number of permutation of i elements, i!. Then

B 2 d [i] = i-1 j=0 d-2j 2 i! . (27)
B.2 Computation of KL (N (0, Σ 1 ) N (0, Σ 2 ))

Let us prove that if Σ 1 and Σ 2 are two covariance matrix, then KL (N (0, Σ 1 ) N (0, Σ 2 )) =

v∈Sp(A) v -log(1 + v) 2 (28)
where

A = (Σ -1 2 -Σ 1 -1 )Σ 1 .
First of all, for a covariance matrix Σ, the density f Σ of N (0, Σ) is given by

∀x ∈ R d f Σ (x) = 1 (2π) d/2 (det Σ) 1/2 exp - 1 2 x T Σ -1 x . ( 29
)
We compute the KL divergence between f Σ 1 and f Σ 2

KL (f Σ 1 f Σ 2 ) = log f Σ 1 (x) f Σ 2 (x) f Σ 1 (x)dx (30) = 1 2 log det Σ 2 det Σ 1 f Σ 1 (x)dx =1 (31) + 1 2 x T Σ -1 2 xf Σ 1 (x)dx =Tr(Σ -1 2 Σ 1 ) (32)
+ 1 2 x T Σ 1 -1 xf Σ 1 (x)dx =Tr(Σ 1 -1 Σ 1 )=d (33) = 1 2 log det Σ 2 -log det Σ 1 + Tr Σ -1 2 Σ 1 -d (34)
We remark that

Tr Σ -1 2 Σ 1 -d = Tr(A) = v∈Sp(A) v (35)

  d d and Part d = Part d d .

2 d

 2 computed in ISDE as estimators of the quantities log(f S )f S S∈Set 2 d

Figure 1 :

 1 Figure 1: Comparison of empirical log-likelihood on validation data for different density estimation methods

Figure 2 :

 2 Figure 2: Comparison of empirical log-likelihood on validation data for best, worst and random partitions

Firstly, we

  can visualize the edit distance from P to the 10 best partitions (excluding P) in the sense of the problem of partition selection, 10 random partitions, and the 10 worst partitions.

Figure 3 :

 3 Figure 3: Edit distance from P for 10 best, 10 random and 10 worst partitions

Figure 4 :

 4 Figure 4: Mean log-likelihood on validation data with respect to edit distance from P for the partitions visited by the random walk

Figure 5 :

 5 Figure 5: Distribution of ( n (P)) P∈Part 5 13

Figure 6 :

 6 Figure 6: Mean log-likelihood on validation data with respect to edit distance from P for 1, 941 best partitions

  implementation 0.0006 0.0020 0.0073 0.0176 0.0684 0.1163 Scikit-learn implementation 0.0008 0.0126 0.1952 1.1564 4.915121.3631

B 2 d

 2 Let us prove the following formula : integer i, let us denote by B 2 d [i] the number of partitions of Part k d with exactly i blocks of size 2. A first remark is that B 2 d [i] = 0 as soon as i > d/2 , then evaluate B 2 d [i]. It is not hard to count the number of possibilities to select i pairs of distinct elements of {1, . . . , d} taking into account in which order there were selected. For the first pair, there are d 2 choices, then d-2 2 choices for selecting another pair among the other variables, and so on. Then there are i-1 j=0 d-2j 2 ordered pairs of variables of {1, . . . , d}.

Table 1 :

 1 S∈ P is a list of bandwidths. For S ∈ Set k d , fS,h S denotes an estimator of the form 2 constructed with the features belonging to S. Let us start by comparing S d and B d , the respective cardinals of Set d and Part d . We have S d = 2 d -1 and B d is known as the Bell number of order d. table1shows how these quantities compare for dimension lying between 10 and 15. Number of partitions vs number of subsets

	Number of Partitions vs. Number of Subsets Before starting the explana-
	tion of how ISDE works, let us highlights some comparison between the number
	of partitions in Part k d and the number of subsets in Set k d .
	d	10	11	12	13	14	15
	S d	1, 023	2, 047	4, 095	8, 191	16, 383	32, 767
	B d 115, 975 678, 570 4, 213, 597 27, 644, 437 190, 899, 322 1, 382, 958, 545
	We remark that the number of partitions is much higher than the number of
	features. Even if we restrict ourselves to small values of k, the difference remains
	important. We denote S k d and B k d the cardinals of Set k d and Part k d . It is simple to
	see that

Table 2 :

 2 × 10 10 6.1 × 10 17 7.3 × 10 25 2.8 × 10 34 Number of partitions vs number of subsets

	1)
	2 ( d/2 )!	(10)

  that the logarithm changes products into sums and the linearity of the operator P [.]. By opposition of what we have seen for the squared L 2 loss, if S ∈ P 1 and S ∈ P 2 , we will have h P 1 S = h P 2 S . Then under KL loss, bandwidths optimization over marginal estimators and partition selection can be decoupled, leading to the necessity of computing S k d density estimators instead of B k d . As shown in table 1 and table 2, it leads to an appreciable gain in terms of algorithmic complexity.

	h P S ) S∈P	P log fP,h P = max (h P S ) S∈P S∈P	P log fS,h S	(14)
	using the property

  It is versatile: even if we present the construction of ISDE using KDEs for marginal estimation, it is possible to use any other base multivariate density estimator.

	Conclusion The resulting algorithm is algorithm 1. It enjoys the following prop-
	erties:
	• It exploits the decoupling of marginal density estimation and partition se-
	lection offered by choice of KL as discrepancy measure: it optimizes over
	partitions in Part k d even if it only requires the computation of Set k d marginal
	estimators.
	d • P∈Part k		S∈P	n (S) using linear programming
	formulation
	end
	Algorithm 1: ISDE
				but less than 10 seconds with LP
	formulation.
	d	9	10 11 12 13	14	15
	Brute-Force Approach 0.2 0.9 5.2 32 219 1304 10437
	LP Solver	0.1 0.2 0.4 0.8 1.9 4.1	9.1
	Table 3: Running time (seconds): linear programming vs brute-force approach for
	partition selection

input : X 1 , . . . , X N ∈ R d , k integer with k ≤ d, integers m and n and a subroutine to perform multidimensional density estimation output: Partition P ∈ Part k d , marginal estimates ( fS ) S∈ P begin for S ∈ Set k d do Compute fS (W 1 , . . . , W m ) thanks to the density estimation subroutine Compute n (S) end Compute P ∈ arg max

Table 4 :

 4 Empirical log-likelihood on validation data for different density estimatorsConclusion For [2, 2, 1], ISDE and FDE give similar results as they output the same graph and the same bandwidths. They both outperform CVKDE.

	ISDE	1.83 ± 0.08 4.05 ± 0.15 6.30 ± 0.25
	FDE	1.83 ± 0.08 2.88 ± 0.14 5.89 ± 0.33
	CVKDE 0.56 ± 0.03 3.49 ± 0.11 3.96 ± 0.16

Table 5 :

 5 Comparison of running time (seconds) of sklearn implementation and ours for KDE constructed on n points and evaluated on n points The computation of the quantities ( n (S)) S∈Set k d requires many repetitions of KDE evaluation. In table 6 we provide estimation of the running time for this step for various values of k and d and considering a 5-fold cross-validation to estimate each bandwidth among 30 candidate values. The quantities m and n are both set to 1, 000.

	k	d 5	10	20	30	40	50
	2	2.0 6.8	18	40	69.15	108
	3	3.3 23	121	409	949	1, 862
	4	3.9 53	590	3, 053	9, 572	23, 536
	5	4.2 61 2, 154 17, 589 75, 228 233, 323

Table 6 :

 6 Running time (seconds) for ( n (S)) S∈Set k d computation with respect to k and d and with 5-fold cross selected bandwidths over 30 possible values and for m = n = 1000 Once the quantities ( n (S)) S∈Set k d are computed, it remains to perform partition selection. As mentioned before we use the python package Pulp [16]. The running time of this step for different values of k and d are presented in table 7.

	k	d	5	10	20	30	40	50
	2		0.02 0.03 0.08 0.20 0.47 0.84
	3		0.02 0.05 0.41 1.9	6.2	15
	4		0.02 0.09 2.0 14.8 62	190
	5		0.02 0.13 7.3 84.7 482 2, 045

Table 7 :

 7 Running time (seconds) for partition selection step with respect to k and d

Table 8 :

 8 Gaussian: KL Losses (.10 3 ) -σ = 0.7, N = 6000 Recovery We are interested not only in performance, but we also want to find the correct partition in order to get qualitative information about datasets. In table 9 we collect, for the same experiment as above, the rate of recovery of the proper partition. In parentheses is displayed the rate of admissible output partition: a partition is admissible if all the blocks of the original partition are subsets of blocks of this one.

	S	[2, 2]	[4, 4, 1]	[4, 3, 2, 3] [4, 4, 3, 3 ,2 ]
	ISDE	0.60 ± 0.21 1.88 ± 0.52 2.85 ± 0.60	5.30 ± 0.96
	BDCS	0.60 ± 0.21 1.72 ± 0.46 2.63 ± 1.01 4.42 ± 1.80
	Empirical 0.80 ± 0.20 3.62 ± 0.53 6.88 ± 0.84	12.63 ± 0.83

Table 9 :

 9 Gaussian: Recovery -σ = 0.7, N = 6000

Acknowledgement This work was supported by the program Paris Region Ph.D. of DIM Mathinnov and was partly supported by the French ANR Chair in Artificial Intelligence TopAI -ANR-19-CHIA-0001. The author is thankful to Marc Glisse and Pascal Massart for their constructive remarks on this work.

Code availability

The code to reproduce the experiments presented here is available at https://github.com/Louis-Pujol/ISDE-Paper.

Data availability Original datasets were downloaded from the repository presented in [START_REF] Lukas | Comparison of clustering methods for high-dimensional single-cell flow and mass cytometry data[END_REF] and available at the address https://flowrepository.org/id/FR-FCM-ZZPH.

We also remark that log

Combining these results with equation (34) leads to the desired formula.


Ilkka Törmä
email: iatorm@utu.fi

Complexity of Generic Limit Sets of Cellular Automata

Keywords: Cellular automata, Limit set, Generic limit set, Topological dynamics

The generic limit set of a topological dynamical system is the smallest closed subset of the phase space that has a comeager realm of attraction. It intuitively captures the asymptotic dynamics of almost all initial conditions. It was defined by Milnor and studied in the context of cellular automata, whose generic limit sets are subshifts, by Djenaoui and Guillon. In this article we study the structural and computational restrictions that apply to generic limit sets of cellular automata. As our main result, we show that the language of a generic limit set can be at most Σ 0 3 -hard, and lower in various special cases. We also prove a structural restriction on generic limit sets with a global period.

Introduction

One-dimensional cellular automata (CA for short) are discrete dynamical systems that act on the set A Z of bi-infinite sequences by a local rule that is applied synchronously at every coordinate. They can be used to model physical and biological phenomena, as well as massively parallel computation.

The limit set of a topological dynamical system (X, T ) consists of those points that can be seen arbitrarily late in its evolution. Limit sets of cellular automata have been studied by various authors from the computational (e.g. [START_REF] Culik | On the limit sets of cellular automata[END_REF][START_REF] Kari | The nilpotency problem of one-dimensional cellular automata[END_REF][START_REF] Borello | Turing degrees of limit sets of cellular automata[END_REF]) and structural (e.g. [START_REF] Maass | On the sofic limit sets of cellular automata[END_REF][START_REF] Ballier | Limit sets of stable and unstable cellular automata[END_REF]) points of view. In [START_REF] Milnor | On the concept of attractor[END_REF], Milnor defined the likely and generic limit sets of a dynamical system. The likely limit set associated to an invariant probability measure µ on X is the smallest closed subset C ⊂ X such that for µ-almost every x ∈ X, all limit points of (T n (x)) n∈N are in C. The generic variant is a purely topological notion that replaces "for µ-almost every x ∈ X" with "for every x in a comeager subset of X". As far as we know, generic limit sets have been studied relatively little in dynamical systems theory.

In [START_REF] Djenaoui | The generic limit set of cellular automata[END_REF], Djenaoui and Guillon studied the generic limit sets of dynamical systems in general and directed cellular automata (arbitrary paths in spacetime diagrams of CA) in particular. They related dynamical properties of a given CA to the structure of its generic limit set in different directions and its relation to the set of equicontinuity points and the limit set. For example, they proved that the generic limit set of an almost equicontinuous CA is exactly the closure of the asymptotic set of its equicontinuity points. They also provide a combinatorial characterization of the generic limit set of a CA, which allows us to study its descriptional and structural complexity and carry out complex constructions. This point of view was not present in [START_REF] Djenaoui | The generic limit set of cellular automata[END_REF], where relatively simple examples of generic limit sets were provided to highlight the main classification results.

As our main result, we prove that the language of a generic limit set of a CA is always Σ 0 3 , and present an example which is complete for this class. If the generic limit set is minimal, then this bound cannot be attained, since its language is Σ 0 2 . We also prove that the dynamics of the CA on its generic limit set must be nontrivial in complex instances: if the CA is eventually periodic or strictly one-sided on the generic limit set, its language is Σ 0 1 or Π 0 2 , respectively, and if the restriction is a shift map, then the generic limit set is chain-transitive. These restrictions are proved by constructing "semi-blocking words" that restrict the flow of information. Finally, we present a structural restriction for generic limit sets: if they consist of a finite number of two-way chain components for the shift map, then they cannot have a global period.

Definitions

Let X be a topological space. A subset of X is comeager if it contains an intersection of countably many dense open sets.

A dynamical system is a pair (X, f ) where X is a compact metric space and f : X → X is a continuous function. We say (X, f ) has trivial dynamics if f = id X . The limit set of f is Ω f = t∈N f t (X). For x ∈ X, we define ω(x) as the set of limit points of the forward orbit (f t (x)) t∈N , and ω(Y

) = y∈Y ω(y) for Y ⊂ X. The realm (of attraction) of a subset Y ⊂ X is D(Y ) = {x ∈ X | ω(x) ⊂ Y }.
The generic limit set of f , denoted ω(f ), is the intersection of all closed subsets C ⊂ X such that D(C) is comeager; then ω(f ) itself has a comeager realm.

We consider one-dimensional cellular automata over a finite alphabet A. The full shift A Z is a compact metric space with the distance function

d(x, y) = inf{2 -n | x [-n,n] = y [-n,n] }, and the left shift map σ : A Z → A Z , defined by σ(x) i = x i+1 , is a homeomorphism. The cylinder sets [w] i = {x ∈ A Z | x [i,i+|w|) =
w} for w ∈ A * and i ∈ Z form a prebasis for the topology, and the clopen sets, which are the finite unions of cylinders, form a basis. We denote [w] = [w] 0 . A subshift is a closed and σ-invariant set X ⊂ A Z . Every subshift is defined by a set F ⊂ A * of forbidden words as X = A Z \ w∈F i∈Z [w] i , and if F can be chosen finite, then X is a shift of finite type (SFT). The language of X is defined as L(X) = {w ∈ A * | [w] ∩ X = ∅}, and we denote L n (X) = L(X) ∩ A n . The order-n SFT approximation of X is the SFT S n (X) ⊂ A Z defined by the forbidden patterns A n \ L n (X). We say X is transitive if for all u, v ∈ L(X) there exists w ∈ A * with uwv ∈ L(X), and mixing if the length of w can be chosen freely as long as it is large enough (depending on u and v). We say X is chain transitive if each S n (X) is transitive. We say X is minimal if it does not properly contain another subshift; this is equivalent to the condition that for every w ∈ L(X) there exists n ∈ N such that w occurs in each word of L n (X).

A morphism between dynamical systems (X, f ) and (Y, g) is a continuous function

h : X → Y with h • f = g • h. If h is surjective, (Y, g) is a factor of (X, f ). A cellular automaton is a morphism f : (A Z , σ) → (A Z , σ). Equivalently, it is a function given by a local rule F : A 2r+1 → A for some radius r ∈ N as f (x) i = F (x [i-r,i+r]
). The pair (A Z , f ) is a dynamical system. Generic limit sets were defined by Milnor in [START_REF] Milnor | On the concept of attractor[END_REF] for general dynamical systems, and were first considered in the context of cellular automata in [START_REF] Djenaoui | The generic limit set of cellular automata[END_REF].

In this article, a Turing machine consists of a finite state set Q with an initial state q 0 and a final state q f , a tape alphabet Γ that is used on a one-way infinite tape together with a special blank symbol ⊥ / ∈ Γ , and a transition rule δ that allows the machine to move on the tape and modify the tape cells and its internal state based on its current state and the contents of the tape cell it is on. Turing machines can decide any computable language and compute any computable function in the standard way.

We give an overview of the arithmetical hierarchy. A computable predicate over N is Π 0 0 and

Σ 0 0 . If φ is a Π 0 n predicate, then ∃k 1 • • • ∃k m φ is a Σ 0 n+1 formula, and conversely, if φ is Σ 0 n , then ∀k 1 • • • ∀k m φ is Π 0 n+1 .
Subsets of N defined by these formulas are given the same classifications, and we extend them to all sets that are in a computable bijection with N. For these sets, we define ∆ 0 n = Π 0 n ∩ Σ 0 n . The computable sets form ∆ 0 1 and the computably enumerable sets form Σ 0 1 . A subshift is given the same classification as its language.

Auxiliary results

We begin with auxiliary results on generic limit sets of cellular automata that are used in several proofs.

Lemma 1 (Proposition 4.11 in [START_REF] Djenaoui | The generic limit set of cellular automata[END_REF]). Let f be a CA. Then ω(f ) is a nonempty f -invariant subshift.

The following result gives a combinatorial characterization for generic limit sets of cellular automata.

Lemma 2 (Corollary of Remark 4.4 in [START_REF] Djenaoui | The generic limit set of cellular automata[END_REF]). Let f be a CA on A Z . A word s ∈ A * occurs in ω(f ) if and only if there exists a word v ∈ A * and i ∈ Z such that for all u, w ∈ A * there exist infinitely many t ∈ N with f

t ([uvw] i-|u| ) ∩ [s] = ∅.
We say that the word v enables s for f . Lemma 3. Let f be a CA on A Z , let n ∈ N, and let [v] i ⊂ A Z be a cylinder set. Then there exists a cylinder set [w] j ⊂ [v] i and T ∈ N such that for all t ≥ T we have f

t ([w] j ) ⊂ [L n (ω(f ))].
Words w with the above property are called ω(f )-forcing, since they force the word f t (x) [0,n) to be valid in ω(f ) whenever w occurs in x at position j. The result intuitively states that any word can be extended into a ω(f )-forcing word.

Proof. Denote A n \ L n (ω(f )) = {u 1 , . . . , u k }. Since u 1 does not occur in ω(f ), Lemma 2 applied to [v] i implies that there exist words

a 1 , b 1 ∈ A * and T 1 ∈ N such that f t ([a 1 vb 1 ] i-|a1| ) ∩ [u 1 ] = ∅ for all t ≥ T 1 . For u 2 we find words a 2 , b 2 ∈ A * and T 2 ∈ N such that f t ([a 2 a 1 vb 1 b 2 ] i-|a2a1| ) ∩ [u 2 ] = ∅ for all t ≥ T 2 . Continuing like this, we obtain a word w = a k • • • a 1 vb 1 • • • b k , a position j = i -|a k • • • a 1 | and a number T = max(T 1 , . . . , T k ) that have the desired property.
Example 1 (Example 5.12 in [START_REF] Djenaoui | The generic limit set of cellular automata[END_REF]). Consider the minimum CA f :

{0, 1} Z → {0, 1} Z defined by f (x) i = min(x i , x i+1 ). We claim that ω(f ) = { ∞ 0 ∞ }.
Proving this directly from the definition is not difficult, but let us illustrate the use of Lemma 2. First, every word 0 n for n ∈ N is enabled by itself: for all u, v ∈ {0, 1} * and t ∈ N we have

f t ([u0 n v] -|u| ) ⊂ [0 n ]. On the other hand, suppose s ∈ L(ω(f )), so that some cylinder set [w] j enables s. Choose u = v = 0 |j|+1 . Then every x ∈ [uwv] j-|u| satisfies x -|j|-1 = 0, so that f |j|+1+|s| (x) [0,|s|) = 0 |s| .
Since a cell can never change its state from 0 to 1, we have s = 0 |s| . Hence the language of ω(f ) is 0 * , and the claim is proved.

Complexity of generic limit sets

From the combinatorial characterization we can determine the maximal computational complexity of the language of the generic limit set.

Theorem 1. The language of the generic limit set of any CA is Σ 0 3 . For any Σ 0 3 set P , there exists a cellular automaton f such that P is many-one reducible to L(ω(f )).

Proof. The condition given by Lemma 2 is Σ 0 3 . For the second claim, since P is a Σ 0 3 set, there is a computable predicate ψ such that P = {w ∈ A * | ∃m ∀m ∃k ψ(w, m, m , k)}. Let M be a Turing machine with state set Q, initial state q 0 ∈ Q, two final states q 1 f , q 2 f ∈ Q, a read-only tape with alphabet Γ A = A ∪ {#, $} and a read-write tape with some tape alphabet Γ with special symbol 1 ∈ Γ . Both tapes are infinite to the right, and M has only one head that sees the same position of both tapes. When initialized in state q 0 , the machine checks that the read-only tape begins with #w#$ m # for some w ∈ A * and m ≥ 0, and the read-write tape begins with 1 3n+5 ⊥ for some n ≥ 0, halting in state q 1 f if this is not the case. Then it enumerates n pairs (m , k) ∈ N 2 , starting from (0, 0) and moving from (m , k) to (m + 1, 0) if ψ(w, m, m , k) holds, and to (m , k + 1) otherwise. If the process ends with k > 0, then M halts in state q 1 f . Otherwise it writes 0s to the |w| + 2 leftmost cells of the read-write tape, goes to the leftmost cell and halts in state q 2 f . Then w ∈ P if and only if for some m ∈ N, the machine halts in state q 2 f for infinitely many choices of n; denote this condition by M (w, m, n).

Denote Σ M = (Q ∪ {←, →}) × Γ and Σ 0 = {B, E, S 1 , S 2 , S 2 , S 3 , }. We construct a radius-3 CA f on the alphabet Σ = (Σ M ∪ Σ 0 ) × Γ A to whose generic limit set P reduces. We write elements of Σ M × Γ A as triples (q, g, a) ∈ (Q ∪ {←, →}) × Γ × Γ A . The first track of f contains elements of Σ M , which are used to simulate computations of M , and Σ 0 , which perform a geometric process that initialized such simulations. The element B forms a background on which the signals E, S 1 , S 2 , S 2 and S 3 travel. The last track of Σ is never modified by f , and it serves as the read-only tape of M in the simulation. We think of f as a non-uniform CA over Σ M ∪ Σ 0 whose local function at each coordinate i ∈ Z depends on the element s ∈ Γ A at i. The automaton f is defined by the following constraints:

-The signal E always travels to the right at speed 2. For k = 1, 2, 3, as long as the signal S k or S k has Bs to its right, it travels to the right at speed k. The signals E, S 1 and S 2 produce Bs in their wake, while S 3 produces (←, 1)-states and S 2 produces (←, ⊥)-states. -When the signals S 2 and S 1 collide, they produce the four-cell pattern (q 0 , 1)(←, 1)S 3 , where S 3 lies at the point of their collision. When S 3 and S 2 collide, they are replaced by an S 2 .

-In an interval of the form (→, g 0 ) . . . (→, g m-1 )(q, g m )(←, g m+1 ) . . . (←, g m+n ) that is either unbounded or terminated on its right by S 3 or S 2 , f simulates a computation of M using q as the head, the Γ -track as the read-write tape and the Γ A -track as the read-only tape. If q = q i f is a final state, it is replaced by E instead.

-Any pattern not mentioned above produces E-states.

In particular, the signals S 1 and S 2 are never created, so they always originate from the initial configuration. The signal S 3 and all Turing machine heads originate either from the initial configuration or a collision of S 2 and S 1 , and S 2 originates from the initial configuration or a collision of S 3 and S 2 . An E-signal, once created, cannot be destroyed.

For a word w ∈ A + , define ŵ = (q 2 f ← |w|+1 , 0 |w|+2 , #w#). We claim that w ∈ P if and only if ŵ ∈ L(ω(f )). The proof is visualized in Figure 1.

Suppose first that w ∈ P , so that there exists m ∈ N such that M (w, m, n) holds for infinitely many n. We claim that the word w = (E |w|+m+3 , #w#$ m #) enables ŵ, so let u, v ∈ Σ * be arbitrary. We construct a configuration x ∈ [u wv] -|u| , which corresponds to the horizontal line of Figure 1, as follows. To the right of v we put (E, #) ∞ , and to the left of u we put only #-symbols on the second track. Let n > |u| + 4 be such that M (w, m, n) holds. On the first track of x [-n+2,-n+4] , put S 1 S 2 E; on x -2n+2 , put S 2 ; on all remaining cells put B. The E-signals will destroy everything in their path and replace them with B-cells, so we can ignore the contents of the first track of x [-|u|,∞) . The

E S2 S1 S2 S3 u w v B B (←, 1) (←, ⊥) q S 2 E
Fig. 1. Proof of Theorem 1, not drawn to scale. Time increases upward.

S 1 -signal and the leftmost S 2 -signal will collide at coordinate 2 of f n (x) (the white circle in Figure 1), resulting in the pattern (q 0 , 1)(←, 1)S 3 at coordinate -1 and Bs to its left. The simulated computation of M begins at this time step. The resulting S 3 -signal collides with the rightmost S 2 -signal at coordinate 3n + 5 of f 2n+1 (x), producing (←, 1)-states until that point and transforming into a S 2 that produces (←, ⊥)-states. This means M has 1 3n+5 ⊥ on its the readwrite tape at the beginning of the computation, and #w#$ m # on the read-only tape. Hence it eventually writes 0 |w|+2 to the tape and halts in state q 2 f at some time step t ∈ N. Then f t (x) ∈ [ ŵ], and we have showed that w enables ŵ, so ŵ ∈ L(ω(f )) by Lemma 2.

Suppose then ŵ ∈ L(ω(f )), so that ŵ is enabled by some word w ∈ Σ * at coordinate i ∈ Z. We may assume i ≤ 0 and |w | ≥ i + |w| + 2 by extending w if needed. Let k ≥ 0 and choose u = (E, #) k and v = (E, #). In a configuration x that contains uw v, any Turing machine head to the right of u is eventually erased by the E-symbols. Those within w are erased after |w | steps, and those to the right of w are erased before they reach the origin. Thus, if f t (x) ∈ [ ŵ] for some t > |w |, then the q 2 f in this configuration is the head of a Turing machine produced at the origin by a collision of some S 2 -signal and S 1 -signal at some earlier time t < t (again the white circle in Figure 1). After a finite computation, M can halt in state q 2 f only at the left end of the tape, so the collision happens at coordinate 2 and t > |w |. Since the signals S 2 and S 1 cannot be created, they originate at coordinates -2t -2 and -t -2 of x. Since the Turing machine eventually halts in state q 2 f , after being initialized it will read 1 3n+5 ⊥ on its read-write tape and #w#$ m # on the read-only tape for some m, n ∈ N with M (w, m, n). Since the read-only tape cannot be modified by f , w already contains the word #w# on its second track, and since v has # on its second track, w must contain #w#$ m # p for some p ≥ 0. Hence m is independent of k.

The signal S 3 produced at the same collision as q 0 continues to the right at speed 3, producing (←, 1)-states until it is destroyed. To its right we have B-states produced by the initial E-signals in u, followed by those E-signals.

Since the Turing machine reads 1 3n+5 ⊥ on its tape, the S 3 -signal is destroyed after n + 1 steps, at time t + n + 1, either by encountering an invalid pattern or by collision with an E-signal or S 2 -signal. In the first two cases, after the removal of S 3 the segment of (←, 1)-states produced by it is now bordered by an E-state, which is an invalid pattern and results in new E-states by the last rule of f . These E-states will eventually destroy the entire computation segment before the Turing machine can halt. Hence S 3 must collide with an S 2 -signal at coordinate 3n+5 at time t +n+1. This signal originates at position -2t +n+3 in x, which must be to the right of the S 1 -signal at coordinate -t -2 that produces S 3 , since these signals do not collide. Hence n > t -5 > k + |w | -5, so n grows arbitrarily large with k. We have shown w ∈ P .

If we know more about the structure of ω(f ) and the dynamics of f on it, we can improve the computability bound.

Proposition 1. Let f be a CA. If ω(f ) is a minimal subshift, then its language is Σ 0 2 .
Proof. Denote X = ω(f ) and let w ∈ L(X). Since X is minimal, there exists n ∈ N such that w occurs in each word of L n (X). Let [v] j be an X-forcing cylinder set that satisfies f t ([v] j ) ⊂ [L n (X)] for all large enough t, as given by Lemma 3. For these t, the set f t ([v] j ) intersects [w] i for some 0 ≤ i ≤ n-|w|. On the other hand, if w / ∈ L(X) then such a word v does not exist, since each word can be extended into one that eventually forbids w. This means that w ∈ L(X) is equivalent to the Σ 0 2 condition that there exist v ∈ A * , j ∈ Z, n ∈ N and T ∈ N such that for all t ≥ T we have f

t ([v] j ) ∩ n-|w| i=0 [w] i = ∅.
Proposition 2. Let f be a CA and suppose that its restriction to ω(f ) is equicontinuous. Then ω(f ) has a Σ 0 1 language.

Proof. Denote X = ω(f ). An equicontinuous CA on any subshift is eventually periodic (this was shown in [START_REF] Kůrka | Languages, equicontinuity and attractors in cellular automata[END_REF] for the full shift, and the general case is not much more difficult), so that there exist k ≥ 0, p ≥ 1 with

f k+p | X = f k | X .
Let r ∈ N be a common radius of f and f p , and let [w] j ⊂ A Z and T ∈ N be given by Lemma 3, so that f t ([w] j ) ⊂ [L 3r (X)] for all t ≥ T . By extending w if necessary, we may assume j ≤ 0, |w| = 3r + 2h and f t (x) [r,2r) = f t (x ) [r,2r) for all x, x ∈ [w] j and t ∈ [T, T + k + p), where h = |j|. Since f has radius r and is eventually periodic on X, we then have

v t := f t (x) [r,2r) = f t (x ) [r,2r
) for all t ≥ T , and the sequence of words (v t ) t≥T +k is p-periodic.

Let n ≥ 0 and u ∈ A 2n . For t ≥ T we have

f t ([wuw] -|w|-n ) ⊂ [v t ] -2r-h-n ∩ [v t ]
n+h+r , so that no information can be transmitted over the v t -words. For x ∈ [wuw] -|w|-n the sequence of words s = (f t (x) [-2r-h-n,n+h+2r) ) t≥T only depends on its values at t ∈ [T, T + k + p), and eventually contains only words of L(X) since we may extend the central pattern of x into one that is X-forcing. Thus s is eventually p-periodic. Since each word s t+1 is determined by s t using the local rule of f , the eventually periodic part is reached when a repetition occurs. The prefixes and suffixes of length r of each word s t already form pperiodic sequences (since they are equal to v t ), so this happens after at most p|A| 2(r+h+n) steps.

Let v ∈ A * be arbitrary. By Lemma 2, v ∈ L(X) if and only if there is a cylinder set

[v ] i with f t ([u v w ] i-|u | ) ∩ [v] = ∅
for all u , w ∈ A * and infinitely many t. By extending v if necessary, we may assume [v ] i = [wuw] -|w|-n for some n ≥ |v| and u ∈ A 2n . Then v occurs infinitely often in words of the eventually periodic sequence s. We have shown that v ∈ L(X) if and only if there exist n ≥ |v| and u ∈ A 2n with f

t ([wuw] -|w|-n ) ∩ [v] = ∅ for some t ∈ [T, T + p|A| 2(r+h+n) + k + p)
, where w, T and h are fixed. Hence X has a Σ 0 1 language.

Proposition 3. If a CA f is the identity on ω(f ), then ω(f ) is a mixing subshift.

Proof. Let v 1 , v 2 ∈ L(X) be arbitrary, and let [w 1 ] i1 , [w 2 ] i2 be two cylinder sets that enable them and are X-forcing in the sense that f t ([w j ] ij ) ⊂ [L n+2r ] -r for all t ≥ T . We may assume, by extending the w j and increasing T if necessary, that f T (x) [0,n) = v j for all x ∈ [w j ] ij , and then f t ([w j ] ij ) ⊂ [v j ] for all t ≥ T . For all large enough N the intersection [w 1 ] j1 ∩ [w 2 ] j2+N is nonempty, and hence contains an X-forcing cylinder

[u] k with f t ([u] k ) ⊂ [L N +|v2| (X)
] for all large enough t. This implies that v 1 A N -|v2| v 2 intersects L(X) for all large enough N , i.e. X is mixing.

We say that a CA f : X → X on a subshift X ⊂ A Z is eventually oblique if f n has a neighborhood that is contained in (-∞, -1] or [1, ∞) for some (equivalently, all large enough) n ∈ N. All shift maps except the identity are eventually oblique.

Proposition 4. Let f be a CA and suppose that its restriction to ω(f ) is eventually oblique. Then ω(f ) has a Π 0 2 language.

Proof. Denote X = ω(f ) and let w ∈ A * . We claim that w ∈ L(X) if and only if the empty word enables w for f , which is a Π 0 2 condition. By Lemma 2 it suffices to prove the forward direction, and the idea of the proof is the following. Since w occurs in X, it is enabled by some cylinder set, which we can extend into one that eventually forces a long segment to contain patterns of X. On this segment, information can flow only from right to left under iteration of f . We use another X-forcing cylinder to block all information flow from the enabling word to the right-hand side of this segment before it is formed. Then the contents of the segment are independent of the word that originally enabled w, so we can swap it for any other word. The argument is visualized in Figure 2.

We assume without loss of generality that some f n has [1, r] as a neighborhood on X, where r ∈ N is a common radius for each f h on A Z for 0 ≤ h ≤ n.

Lemma 3 gives us a cylinder set [v]

i ⊂ A Z and T ∈ N with f t ([v] i ) ⊂ [L 3r (X)]
for all t ≥ nT . By extending v if necessary, we may assume -|v| ≤ i ≤ -rnT .

Assume w ∈ L(X), so that there exists a cylinder set [u] j ⊂ A Z that enables it for f . For the empty word to enable w, it suffices to show that for an arbitrary cylinder set [u ] j there exist infinitely many t ∈ N with f t ([u ] j ) ∩ [w] = ∅. By extending u and/or u if necessary, we may assume j = j and |u| = |u | ≥ |j|, and denote k = j + |u|. Consider the cylinder set [uv] j . By Lemma 3, there exists a cylinder set [auvb] j-|a| with f t ([auvb] j-|a| ) ⊂ [L 3r+k-i (X)] -2r for all large enough t ∈ N. For the same reason, there exists another cylinder set [cau vbd] j-|ca| with f t ([cau vbd] j-|ca| ) ⊂ [L 3r+k-i (X)] -2r for all large enough t ∈ N. Let T ≥ T be such that nT is a common bound for these conditions. Denote K = nT + n|uvbd| -j, and let t ≥ K be such that f t (x) ∈ [w] for some x ∈ [cauvbd] j-|ca| . There are infinitely many such t since u enables w for f . Let y ∈ [cau vbd] j-|ca| be the configuration obtained by replacing the u in x by u .

We claim

f t (y) ∈ [w]. Note first that f hn (x) [k+sr,∞) = f hn (y) [k+sr,∞)
for all h ≤ T since this holds for h = 0 and r is a radius for f n . This is represented in Figure 2 by the lower dashed line. Since x, y ∈

[v] k , for all s ≥ nT we have f s (x), f s (y) ∈ [L 3r (X)] k-i . From i ≤ -rnT it follows that k -i ≥ k + rnT , so that f nT (x) and f nT (y) agree on [k -i, ∞). Since f n has [-r,
r] as a neighborhood on A Z and [1, r] as a neighborhood on X, for each j ∈ [k -i+r, ki + 2r] the value of f s+n (z) for z ∈ {x, y} depends only on f s (z) [ +1, +r] . Thus if f s (x) and f s (y) agree on [k -i + r, ∞) for some s ≥ nT , then so do f s+n (x) and f s+n (y). By induction, we obtain f ns (x) [k-i+r,∞) = f ns (y) [k-i+r,∞) for all s ≥ T . Denote g(s) = max(-r, k -i + r -s). Then f n(T +s) (x) and f n(T +s) (y) agree on [g(s), ∞) for all s ≥ 0. This is represented in Figure 2 by the upper dashed line. For s = 0 this is true by the previous paragraph, so suppose it holds for some s ≥ 0. Since x ∈ [auvb] j-|a| and y ∈ [cau vbd] j-|ca| , we have f n(T +s) (x), f n(T +s) (y) ∈ [L 3r+k-i (X)] -2r . As in the previous paragraph, the value of f n(T +s+1) (z) for ∈ [g(s) -1, k -i] and z ∈ {x, y} depends only on f n(T +s) (z) [ +1, +r] . The claim follows by induction.

Writing t = ns + h for 0 ≤ h < n, the configurations f t (x) and f t (y) agree on [g(s) + r, ∞) = [0, ∞) since r is a radius for f h . In particular

f t (y) ∈ [w], so that f t ([u ] j ) ∩ [w] = ∅.
Proposition 5. If the restriction of a CA f to ω(f ) is a shift map, then ω(f ) is a chain transitive subshift.

Proof. Suppose f | X = σ n | X for some n ∈ Z. By symmetry we may assume n > 0. Let v 1 , v 2 ∈ L(X) be two words of equal length m, and let N ≥ m be arbitrary. We claim that there exist M > 0 and words u 0 , . . . , u M ∈ L N +n (X) such that v 1 is a prefix of u 0 , v 2 is a prefix of u M and the length-N suffix of each u i is a prefix of u i+1 , which implies the chain transitivity of X. For this, let [w] j be an X-forcing cylinder with f t ([w] j ) ⊂ [L 2r+N +n (X)] -r for all large enough t. By the proof of Proposition 4, [w] j enables both v 1 and v 2 , so there exist x ∈ [w] j and T ≤ t 1 < t 2 with f ti (x) [0,m) = v i for i = 1, 2. Choose M = t 2 -t 1 and u k = f t1+k (x) [0,N +n) . Since r is a radius for f , these words have the required properties.

Using these results, we can prove that some individual subshifts cannot occur as generic limit sets.

Example 2. There is no CA f : A Z → A Z whose generic limit set is the orbit closure of ∞ 01 ∞ . Suppose for a contradiction that there is one. Since X = ω(f ) is invariant under f , we have f ( ∞ 01 ∞ ) = σ n ( ∞ 01 ∞ ) for some n ∈ Z, and then f | X = σ n | X . If n = 0, then Proposition 3 implies that X is mixing, and if n = 0, then Proposition 5 implies that X is chain transitive, but it is neither.

By the results of [START_REF] Olli | Endomorphisms of Sturmian systems and the discrete chair substitution tiling system[END_REF], all cellular automata on Sturmian shifts are restrictions of shift maps. Propositions 2 and 4 imply that the language of a Sturmian generic limit set is Π 0 2 , and Proposition 1 (or the folklore result that every minimal Π 0

n subshift is Σ 0 n ) implies that it is Σ 0 2 .
Hence we obtain the following. Corollary 1. If a Sturmian shift is the generic limit set of a CA, then its language is ∆ 0 2 .

Periodic factors

In some situations, a nontrivial finite factor forbids a subshift from being realized as a generic limit set.

Definition 1. Let X ⊂ A Z be a subshift. The chain relation of width n is the relation on L n (X) defined by u ∼ n v if there exists x ∈ X with x [0,n) = u and x [k,k+n) = v for some k ≥ 0. The symmetric and transitive closure of ∼ n is the σ ± -chain relation of width n. If each ∼ n is equal to L n (X) 2 , we say X is σ ±chain transitive. A σ ± -chain component of X is a maximal σ ± -chain transitive subshift of X.

It is not hard to see that every subshift is the union of its σ ± -chain components, which are disjoint. SFTs and sofic shifts have a finite number of such components, but in other cases their number may be infinite.

Example 3. Let X ⊂ {0, 1, 2} Z be the union of the orbit closures of ∞ 02 ∞ and ∞ 12 ∞ . For each n ∈ N, we have 0 n 2 n , 1 n 2 n ∈ L(X), which implies 0 p 2 n-p ∼ n 2 n and 1 p 2 n-p ∼ n 2 n for all 0 ≤ p ≤ n. Since L n (X) consists of exactly these words, X is σ ± -chain transitive. Lemma 4. Let f be a CA on A Z such that ω(f ) has a finite number of σ ±1chain components X 1 , . . . , X k . Then there is a cyclic permutation ρ of {1, . . . , k} such that f (X i ) = X ρ(i) for each i = 1, . . . , k.

Proof. Since the image of a σ ±1 -chain transitive subshift by a cellular automaton is also σ ±1 -chain transitive, each X i is mapped into some other component X ρ(i) . This defines a function ρ : {1, . . . , k} → {1, . . . , k}. By Lemma 1, ρ is surjective, hence a permutation, and f maps each X i surjectively to X ρ(i) . Corollary 4.13 in [START_REF] Djenaoui | The generic limit set of cellular automata[END_REF] implies that ρ must be a cyclic permutation. Proposition 6. If a subshift X ⊂ A Z has a finite number of σ ±1 -chain components and a finite factor that does not consist of fixed points, then it is not the generic limit set of any CA.

Proof. Suppose that f is a CA with ω(f ) = X, and let π : (X, σ) → (Y, g) be the morphism onto a nontrivial finite factor. Let X 1 , . . . , X k be the σ ±1 -chain components of X. By Lemma 4 we may assume f (X i ) = X i+1 mod k for all i. If Y p ⊂ Y is the subsystem of p-periodic points, then π -1 (Y p ) ⊂ X is an f -invariant subshift consisting of σ ±1 -chain components, so it is nonempty for exactly one p, and p > 1 by the assumption that Y does not consist of fixed points. By taking a factor map from Y onto Z p if necessary, we may assume (Y, g) = (Z p , +1) where the addition is modulo p. Define q : X → Z p by q(x) = π(x) -π(f (x)). Then q is continuous and shift-invariant, and is constant in each component X i . Let r be a common radius for f and right radius for π and q, meaning that π(x) and q(x) are determined by x [0,r) . The right radii exist since π has some two-sided radius s by continuity and satisfies π(x) = π(σ s (x)) -s, which is determined by x [0,2s] , and similarly for q. For w ∈ L 2r+1 (X), denote π(w) = π(x) and q(w) = q(x) for any x ∈ [w] -r .

Let u ∈ L 3r (X 1 ) be arbitrary, let [v] i be a cylinder that enables it given by Lemma 2, and let [w] j ⊂ [v] i be an X-forcing cylinder given by Lemma 3, so that f t ([w] j ) ⊂ [L 3r (X)] for all t ≥ T . We may assume that f T ([w] j ) ⊂ [u] by extending w if necessary. Let m ∈ N be such that [w] j ∩ [w] j+pm+1 is nonempty, which holds for all large enough m. Finally, let [v ] ⊂ [w] j ∩ [w] j+pm+1 be an X-forcing cylinder with f t ([v ] ) ⊂ [L 3r+pm+1 (X)] for all large enough t.

Let x ∈ [v ] be arbitrary. Then we have f T (x) ∈ [u] ∩ [u] pm+1 and f t (x) ∈ [L 3r (X)]∩[L 3r (X)] pm+1 for all t ≥ T . This implies π(f t+1 (x) [r,2r) ) = π(f t (x) [r,2r) )+ q(f t (x) [r,2r) ) and π(f t+1 (x) pm+1+[r,2r) ) = π(f t (x) pm+1+[r,2r) )+q(f t (x) pm+1+[r,2r) ) for each t ≥ T . Since u ∈ L 3r (X 1 ), f permutes the components X i and q is constant in each component, we have π(f t (x) [r,2r) ) = π(f t (x) pm+1+[r,2r) ) for t ≥ T . For large enough t, this is a contradiction with the fact that w (t) = f t (x) [0,pm+3r] ∈ L(X) satisfies π(w As a corollary, a transitive but nonmixing SFT cannot be the generic limit set of a CA, since it admits a factor map to some Z p with p > 1.

Future work

In this article we presented a construction for a generic limit set with a maximally complex language. We also showed that structural and dynamical properties may constrain this complexity, but did not prove the strictness of these bounds. Furthermore, we showed that in some cases a global period is an obstruction for being a generic limit set, but did not prove whether this always holds. Using a construction technique presented in [START_REF] Boyer | µ-limit sets of cellular automata from a computational complexity perspective[END_REF] to study µ-limit sets, we believe we can answer these and other questions.

A cellular automaton f is nilpotent, meaning that f t (A Z ) is a singleton set for some t ∈ N, if and only if its limit set is a singleton [START_REF] Culik | On the limit sets of cellular automata[END_REF]. As a variant of nilpotency, one may define f to be generically nilpotent if ω(f ) is a singleton, as in Example 1. We believe the techniques of [START_REF] Boyer | µ-limit sets of cellular automata from a computational complexity perspective[END_REF] can also be used to characterize the complexity of deciding generic nilpotency.

Fig. 2 .

 2 Fig.2. Proof of Proposition 4, not drawn to scale. Time increases upward. In the shaded region, information flows only from right to left. The configurations f s (x) and f s (y) agree on the part that is right of the dashed line.

  r) ) + i for all i ∈ [0, pm + 2r].