Ément Cl

Richefort

Matthieu Lecouvez

Rob Falgout

Pierre Ramet

TOWARD A MULTILEVEL METHOD FOR THE HELMHOLTZ EQUATION *

Keywords: Multigrid, Helmholtz, Linear Algebra, Polynomial Smoother

It is well known that multigrid methods are very competitive in solving a wide range of SPD problems. However achieving such performance for non-SPD matrices remains an open problem. In particular, two main issues may arise when solving a Helmholtz problem. Some eigenvalues become negative or even complex, requiring the choice of an adapted smoothing method for capturing them. Moreover, since the near-kernel space is oscillatory, the geometric smoothness assumption cannot be used to build efficient interpolation rules. We present some investigations about designing a method that converges in a constant number of iterations with respect to the wavenumber. The method builds on an ideal reduction-based framework and related theory for SPD matrices to correct an initial least squares minimization coarse selection operator formed from a set of smoothed random vectors. We also present numerical results at the end of the paper.

1. Introduction. The numerical simulation of various physical phenomena gives rise to potentially very large linear systems of equations written Au = b in matrix form. These systems can be solved directly by a convenient factorization of A, or iteratively by computing and refining an approximation of the solution u starting from an initial guess u 0 . Multigrid methods [START_REF] Strang | Multigrid methods[END_REF][START_REF] Falgout | An introduction to algebraic multigrid[END_REF] work iteratively and are known to be scalable and quasi-optimal for solving sparse linear systems of equations for many classes of problems. To simplify the discussion in what follows, we use the term "small/large eigenvector" to mean an eigenvector with small/large eigenvalue. We similarly say "positive/negative eigenvector" when referring to the eigenvalue sign.

1.1. General aspects on multigrid methods. The core idea in multigrid methods is to accelerate the computation of u by way of a hierarchy of coarse problems A l u l = b l , l being the level in the grid hierarchy. A restriction operator R l transfers the information from a level l to a coarser one l + 1, while an interpolation operator P l transfers the information from level l + 1 to l. In most symmetric applications, R l = P T l and coarse matrices are constructed following the Galerkin formula A l+1 = P T l A l P l . Two-grid methods actually need both types of solvers: a direct method for the coarse correction, and an iterative method called a smoother on the fine level.

The error propagation matrix for the coarse correction is

(1.1) E = I -P (P T AP ) -1 P T A.

The error propagation matrix for the smoother is

(1.2) E M = I -M -1 A
where M -1 is an approximation of A -1 , for instance the diagonal inverse or a polynomial of A. The smoother is applied before each restriction and after each interpolation. Finding a smoother and a coarse correction that are complementary is a major concern in the design of the method. The interpolator must propagate the coarsest information back to the finest, and transferred errors should be eliminated by the smoother. The smoother targets the large eigenvectors and the coarse correction targets the small eigenvectors. The near-kernel space of smallest eigenvectors is especially important in the design of interpolation. A multi-level method can be created by recursively applying the two-level method to solve the coarse system, where a direct solver is used on the coarsest level. The context in which a multigrid method is applied determines what kind of operators should be used in the method. In elliptic problems, where the convergence of multigrid methods is well known, the matrix A is symmetric positive definite, so smoothers like weighted-Jacobi or Gauss-Seidel are known to be good smoothers since they damp the high frequencies without modifying the low frequency eigenvectors. Likewise, interpolators are designed to target slowly varying components as in classic algebraic multigrid methods (AMG) [START_REF] Stüben | Algebraic multigrid (amg). an introduction with applications[END_REF].

1.2. Why Helmholtz problems are difficult for multigrid. The Helmholtz problem involves indefinite matrices with potentially wide and oscillatory near-kernel space [START_REF] Ernst | Why it is difficult to solve helmholtz problems with classical iterative methods[END_REF]. This complication breaks the near-kernel space geometric smoothness assumption, a keystone of many multigrid methods. To satisfy the complementarity principle in this context, interpolation rules must reproduce the near-kernel oscillation, and smoothers have to deal with both positive and negative eigenvalues. More importantly, finding a recurring process to build a scalable multilevel method is still an open question. The Helmholtz equation (1.3) is our target in this paper.

(1.3) (Continuous Helmholtz problem) ⇔ -∆u -k 2 u = f on Ω + b. c. on ∂Ω
Since the Helmholtz equation can be seen as a shifted Poisson equation, the geometrically smooth components (ie. low Fourier modes) become negative. Because of the shift, the smallest eigenvectors are higher in frequency. Multigrid interpolators must now focus on this more oscillatory spectrum interval. Multiple correction [START_REF] Livshits | Multiple galerkin adaptive algebraic multigrid algorithm for the helmholtz equations[END_REF] and wave-ray [START_REF] Brandt | Wave-ray multigrid method for standing wave equations[END_REF][START_REF] Livshits | A scalable multigrid method for solving indefinite helmholtz equations with constant wave numbers[END_REF] approaches have already been investigated to address this issue. In this paper, we present an approach built on ideal reduction-based ideas, and demonstrate its potential for solving the Helmholtz problem in constant iteration count independent of the wavenumber k. In Section

q(A H A)v = q(σ 2 )v = (1 -p(σ 2 )σ 2 )v
According to Equation (2.2), the lower the amplitude of a component, the less the polynomial smoother will damp it. This intact portion of low components enables the construction of a set of smoothed vectors approximating the near-kernel space, as detailed in Section 3. Let L be a set of points lying in this interval such that

(2.3) ∀ x i ∈ L = (x 1 , . . . , x d ) , x i ̸ = 0 , q(x i ) = 0 ⇔ p(x i ) = 1 x i .
Then the polynomial smoother p is constructed following the Lagrangian formula (2.4)

x i ̸ = x j , p(x) = d j=0 1 x j d i=0,i̸ =j x -x i x j -x i .
However, such interpolation points should not be selected randomly within the interval, but in order to minimize the polynomial amplitude and avoid Runge's phenomenon. Those best interpolating points are defined by the scaled first kind Chebyshev polynomial roots (2.5)

x i := x max + x min 2 + x max -x min 2 cos (2i -1)π 2d .
Large intervals I require a higher polynomial degree to flatten its oscillations.

Constructing an appropriate target interval I. One way to determine

a good interval without preliminary information [START_REF] Baker | Multigrid smoothers for ultraparallel computing[END_REF][START_REF] Adams | Parallel multigrid smoothing: polynomial versus gauss-seidel[END_REF] is to compute a few power iterations to determine x max , overestimate the result by 10%, and choose the lower bound

x min according to x max , for example x min = 0.5x max . To respect the complementarity principle, the percentage of damped eigenvalues by the smoother must approximate the proportion of non-coarse variables. For instance, if a coarse level space is defined by a selection of a quarter of the finer level's degrees of freedom, then three-quarters of the largest amplitude components should be damped by smoothing steps, while the coarse correction deals with the other part. Consequently, since eigenvalues are not necessarily uniformly separated, x min should be determined so that such a proportion of eigenvalues belongs to the interval I. In this paper and following [START_REF] Lin | Approximating spectral densities of large matrices[END_REF], we

This manuscript is for review purposes only.

first compute a rough approximation of the matrix spectral density defined by the distribution function ϕ(t), which represents the probability of finding an eigenvalue at each point t of a given interval. We set the lower bound x min of the Chebyshev nodes interval in a second step so that the probability within the interval is equal to the target proportion, for instance half of the total area in a scenario of exact balance between C and F points. The squared singular values should be scaled by setting

B := 2 xmax A H A -I,
ϕ(t) = ∞ k=0 µ k T k (t) ≈ M k=0 µ k T k (t) , µ k = 2 -δ k0 nπ × Trace(T k (B))
where T k (t) = cos(k arccos(t)) and the coefficients µ k are determined by a moments matching procedure. Here, n corresponds to the matrix size and δ k0 the Kronecker symbol. The trace can be estimated using a set U of n vec random and orthogonal vectors u, where each element of these vectors is chosen following a normal distribution with zero mean and a unit standard deviation. Any vector u can be written as a linear combination of B eigenvectors v, giving the expression u

= n j=1 β j v j . As a consequence, E[u] = 0 and E[β i β j ] = δ ij , giving the following convenient property (2.7) E[u T T k (B)u] = E[ n j=1 β 2 j T k (λ j (B))] = n j=1 T k (λ j (B)) = Trace(T k (B)).
According to (2.7), each trace can be estimated by a sample mean of n vec products

u T T k (B)u,
3 -2 -1 0 1 2 3 4 5 6 0 0.5 1 x 50% x min x max x q(x) = 1 -p(xx)xx q x i Fig. 2.

Spectrum of the polynomial smoother error propagation matrix

This manuscript is for review purposes only.

3. Constructing good interpolation rules. Interpolators are used both to construct the coarse level matrices and to transfer information across levels. SPD and geometric smoothness assumptions cannot be used to determine appropriate interpolators in our case. Some methods such as smoothed aggregation [START_REF] Ek | Convergence of algebraic multigrid based on smoothed aggregation[END_REF][START_REF] Olson | Smoothed aggregation for helmholtz problems[END_REF] and bootstrap-AMG [START_REF] Brandt | Bootstrap amg[END_REF] use candidate vectors that are close to the near-kernel space to design the interpolation rules. These test vectors are either deduced from geometric information [START_REF] Brandt | Wave-ray multigrid method for standing wave equations[END_REF] or algebraically as in adaptive multigrid methods [START_REF] Brezina | Adaptive smoothed aggregation (αsa)[END_REF]. Here, we prefer to stick to a fully algebraic and recurring process to create our interpolators. Candidate vectors will be generated from random vectors smoothed by the polynomial presented in Section 2, and used by the least squares minimization framework to determine good fine variable interpolation rules. This initial least squares interpolator is used as a coarse selection operator in the ideal reduction-based framework [START_REF] Falgout | On generalizing the amg framework[END_REF].

3.1. Ideal framework. Even though the ideal framework in [START_REF] Falgout | On generalizing the amg framework[END_REF] requires an SPD assumption and has not been generalized to indefinite problems, removing orthogonal information from the interpolation range still improves its accuracy. Furthermore, assuming the smoother captures this orthogonal information, we can still expect good convergence as exposed by Equation (3.3) below. Following [START_REF] Falgout | On generalizing the amg framework[END_REF], let C and F be complementary coarse and fine variables subsets of Ω. Let R T : R n C → R n and S : R n F → R n be coarse and fine selection operators respectively, such that RS = 0, for instance, the orthogonal matrices

(3.1) R T = [ 0 I C ] T , S = [ I F 0 ] T .
The space defined by the coarse selection operator R T must be handled by the coarse correction, whereas the fine variables selection operator S defines a space where smoothing must operate in order to respect the complementarity principle. The Ideal Interpolator is a theoretical operator that is the best interpolator satisfying RP = I C in the sense that it minimizes the difference between variables and interpolated coarse variables, within a space that is the most complementary to the range of the smoother.

For this reason, P * is called ideal and is defined by

(3.2) P * = (I -S(S T AS) -1 S T A)R T .
The left operator in (3.2) removes all complementary F-related information from R T . Such information is irrelevant at a coarse level and should be handled by the smoother. Under the assumption that the smoother captures F-related information, the best coarse matrix, according to the complementarity principle, is a matrix where fine variable information is removed. Furthermore, a simple development of a two-grid cycle combining the ideal coarse error propagation operator E * , as initially defined in (1.1), with an F-relaxation matrix E F as a pre-smoother shows that

(3.3) E * E F = (I -P * (P T * AP * ) -1 P T * A)(I -S(S T AS) -1 S T A) = 0.
In (3.3), the final error propagation matrix is null, meaning that in this idealistic scenario, one iteration of F-relaxation and one coarse correction is equivalent to a direct method. As shown in [START_REF] Falgout | On generalizing the amg framework[END_REF], we can extend this principle by separating coarse and fine spaces respectively in directions of low and high frequencies. Let V 0 = (v 1 , . . . , v Card(C) ) and V + = (v Card(C)+1 , . . . , v n ) respectively be low frequency and high frequency eigenvectors sets, assuming here that λ i ≤ λ i+1 . Then we define R T = V 0 and its counterpart S = V + . Naturally, since eigenvectors are orthonormal, the necessary condition RS = 0 is still satisfied. It finally gives This manuscript is for review purposes only.

This example proposes another idealistic dichotomy enabled by P * , maximizing the complementarity principle. The near-kernel space is solved directly at the coarsest level while the high frequencies remain in the smoothing space. Even if (S T AS) -1 is most of the time impossible to use in practice, it gives insight on an idealistic convergence scenario.

Least Squares Minimization

Interpolator. As mentioned at the beginning of Section 3.1, demonstrating that interpolator (3.2) is ideal in the theoretical framework of [START_REF] Falgout | On generalizing the amg framework[END_REF] requires A to be symmetric positive-definite. However, the reduction viewpoint of Equation (3.3) is still valid, hence removing the orthogonal fine information from the coarse selection operator is a viable approach. Numerical experiments show that the coarse selection operator in (3.1) is not a good option for Helmholtz.

Using the lowest components V 0 from Section 3.1 to guarantee the representation of the near-kernel space within the interpolation range is not practical. Instead, we construct a set of vectors approximating an oscillatory and potentially large near-kernel space by using the normal equations polynomial smoother developed in Section 2.

In this section, we present a coarse selection operator R H constructed by a least squares minimization strategy [START_REF] Brandt | Bootstrap amg[END_REF]. Let K be a set of κ vectors that approximate the near-kernel space, and assume some C/F splitting with n C and n F their respective size. Coarse variables are interpolated to the finer level with a simple injection rule, meaning the coarse interpolation block in RH corresponds to a n C × n C identity matrix, while fine interpolation rules are determined by the least squares minimization method presented in this section. Let i be a fine variable and r i the i th row of RH . The idea consists of constructing each fine interpolation rule by minimizing the squared difference between fine values of the near-kernel candidate vectors and the interpolation from their connected coarse variables C i . Denote by K :,l the l th test vector, K i,: a row vector containing the i th values of each test vector, and K Ci,l a vector containing the values in K :,l of the coarse variables that are connected to the i th fine variable. Then

(3.5) ∀i ∈ F , r i = arg min r κ l=1 w l (K i,l -r • K Ci,l ) 2 := arg min r L i (r)
where w l are scaling weights (for instance w l = 1 /λ l if K contains near-kernel eigenvectors). Finding the minimum of the convex loss function L i is equivalent to solving

(3.6) ∇L i (r i ) = 0.
Equation (3.6) can be rewritten element-wise

(3.7) ∂L i (r i ) ∂r ij = κ l=1 2w l (K i,l -r i • K Ci,l )K Cij ,l = 0 , ∀j ∈ [1, card(C i )].
Finally, (3.7) leads to a system of linear equations to solve for each fine variable i

(3.8) r i K Ci W K H Ci = K i W K H Ci
The matrix is full rank and the solution of Equation (3.8) is unique if we have at least κ = max i {Card(C i )} locally linearly independent test vectors. Even if it is statistically always the case when starting from random candidate vectors, the matrix singularity can be detected during the factorization. In that special case, a pseudo-inverse can be computed to find an optimal solution in the least squares sense.

This manuscript is for review purposes only.

(3.9) RH = [ R F I C ] T , Ŝ = [ I F -R H F ] T ,
where RH is the least squares operator presented in Section 3.2 and R F is its fine variable interpolation block. Note that R Ŝ = 0 as required. To simplify the discussion, let A F := ŜH A Ŝ. Beyond the necessity to find appropriate coarse and fine selection operators, another important concern is related to the inverse of A F required in (3.2).

Reorganizing the definition of the ideal interpolator, (

P := (I -ŜA -1 F ŜH A) RH ≈ RH -ŜX -1 K ŜH A RH , 3.10)
where X -1 K is the best polynomial approximating A -1 F within the Krylov subspace K. Since we need to ensure our interpolator keeps good sparsity, the subspace must be constrained according to a given pattern P as detailed in [START_REF] Olson | A general interpolation strategy for algebraic multigrid using energy minimization[END_REF]. In our case, the matrix inversion approximation will be computed column-wise, giving more flexibility than by computing the global constrained matrix approximation at once with a single polynomial. Consequently, let P i be some vector sparsity pattern, and define the associated operator Z i : C n → C Card(Pi) filled with ones and zeros that restricts any full vector to the non-zero pattern P i . Also, define

i := ŜH A RH :,i . (3.11) b
In practice, we choose P i := P(b i ) such that Z T Z i b i = b i . Then, we construct for each right hand side b i the corresponding constrained Krylov subspace (3.12)

K m Pi = Z i b i , Z i A F Z T i Z i b i , . . . , (Z i A F Z T i ) m-1 Z i b i ,
within which we approximate the multiplication of A -1 F with b i . Since each subspace is constructed under a sparsity constraint, the solution is approximated locally with a window of A F . Consequently, it is not guaranteed that the accuracy will increase with respect to m, however a few iterations are enough to reach a good approximation in practice. It is still possible to converge toward the best solution in a least squares sense, but this requires reformulating the problem with normal equations which increases the cost of construction.

κ l = 5 × max i∈F {Card(C i )} + κ l-1 + 10l , κ -1 = 0
This manuscript is for review purposes only.

where C i is the set of coarse variables strongly connected to the fine variable i. The strongly connected variables are selected according to a θ-rule [START_REF] Stüben | Algebraic multigrid (amg). an introduction with applications[END_REF] comparing matrix entries in absolute value. Before being smoothed by the normal equation polynomial, test vectors are created either by the restriction of finer level test vectors or randomly.

In practice, Card(C i ) never exceeds 10. The size m of the Krylov subspaces needed to approximate A -1 F in the ideal framework is set to 3. The three V-cycle schemes benchmarked in Figure 4.1 converge in a roughly constant number of iterations independent of the wavenumber k. We also present these results together with the overall complexity of our multigrid method. Let nnz(•) be the average number of non-zeros per row of a given matrix. We measure the sparsity of each interpolator Pl and level matrix A l in Table 4.1. As expected, matrices on deeper levels are denser. Even if A -1 F is approximated under pattern constraints and already allows to find an interesting trade-off between sparsity and good interpolation properties, it will be necessary to improve the sparsity of P to reach a more competitive multigrid method. A thresholding strategy or other heuristics on its pattern will be one of the main concerns for the future developments of the method. This manuscript is for review purposes only.

Large shift experiment.

Even if the approach presented in this paper can be improved in many ways, it provides a direction for constructing interpolation for problems with oscillatory near-kernel spaces like the Helmholtz equation. However, since the problem is indefinite, the matrix A does not provide a norm. As a consequence, there is no guidance on the convergence. Moreover, even under the convenient assumption where P P T v ≈ v where v is an eigenvector associated with a very small eigenvalue, for some large shift problems, the method can be divergent. To illustrate this issue, let L s be a scaled laplacian matrix shifted by a large coefficient α (yielding λ(L s ) ∈ [-α, 8 -α]), such that its near-kernel space is very oscillatory.

.2 shows the layering of the lowest eigenvector of L s (v -green), the interpolation of the restricted lowest eigenvector ( P P T v -blue), and the vector returned by the coarse correction ( P ( P T A P ) -1 P T Av -red). As we can see, v and P P T are very close in both large shift experiments, the very oscillatory near-kernel vector v is well approximated by the interpolation range. However, in the second experiment, this is not true for the coarse correction where the red vector seems to be oriented oppositely to v. This experiment shows that capturing the oscillatory near-kernel space with an appropriate set of interpolators, such as P , will not necessarily be enough to reach a perfectly recurring method for solving indefinite problems like Helmholtz. For this reason, adding more levels has been challenging. One of our priorities is about finding a framework guaranteeing the convergence of the method in an indefinite context such as this large shift experiment.

This manuscript is for review purposes only.

5.

Conclusions. Indefinite and oscillatory problems are difficult for multigrid methods. The negative eigenvalues require an adapted smoother, and the interpolator should be able to propagate the oscillatory near-kernel space. We presented a method that reaches those requirements up to a certain limit. The normal equation polynomial smoother is designed to target a desired proportion of components according to their amplitude, and the range of our interpolator offers a good approximation of the near-kernel space despite its oscillations. Finding more accurate interpolation rules, improving the sparsity of our operators, and constructing a polynomial without resorting to normal equations will be important points in our future investigations.

However, the ultimate objective is to find a proper framework for indefinite problems guaranteeing the convergence of our multigrid method.

  and the M vectors T k (B)u can be computed from the three-term recurrence relation T k+1 (B)u = 2BT k (B)u -T k-1 (B)u, with initial terms T 0 (B)u = u, T 1 (B)u = Bu. Once the distribution function ϕ is approximated following the expression given in (2.6), a rough area approximation by trapezoid rule yields a correct lower bound x satisfying a proportion around card(F ) /n. Computing x min by remapping x to the initial squared singular values scaling leads to a satisfying and purely algebraic interval in which our polynomial smoother will be the most efficient. The bounds x min and x max are represented in Figure 2.1, where x 50% represents the theoretical lower bound targeted.

-

(3. 4 )

 4 P * = R T and A C = Diag(λ 0 , . . . , λ Card(C) ).

3. 3 .

 3 Ideal approximation from least squares coarse selection. In Section 3.2, we presented a better coarse selection operator for Helmholtz designed by a least squares minimization strategy. Using the framework presented in 3.1, let us define new coarse and fine selection rules

4 .

 4 Numerical Experiments. In this section, we present some numerical results performed with this new multigrid cycle. One pre-smoothing and one post-smoothing iteration of the normal equations Chebyshev polynomial presented in Section 2 are computed before each restriction and after each interpolation respectively. The interpolator P is constructed following Section 3. In practice, the degree d of the normal equation polynomial smoother p of Equation 2.2 is equal to 3, the x max is computed by several power iterations (5 to 10 iterations are enough in practice) and x min by the spectral density approximation described in Section 2.2, with the parameters M = 5 and n vec = 15. The number κ l of input random vectors per level l to construct a correct approximation of near-kernel space follows the arbitrary recursive equation(4.1)

4. 1 .Fig. 4 . 1 .

 141 Fig. 4.1. Number of iterations following the wavenumber k

Fig. 4 . 2 .

 42 Fig. 4.2. Layering of lowest eigenvector (green), interpolation of restricted lowest eigenvector (blue) and coarse correction applied to lowest eigenvector (red) for α = 2.68 (top) vs. α = 2.98 (bottom)

  so the distribution function ϕ can be approximated by a linear

	combination of orthogonal Chebyshev polynomial functions
	(2.6)

Table 4 . 1

 41 Sparsity measurement of level matrices and interpolators following k

	k	70	100	130	160	190	220	250	280	310
	n (×10 5 ) 1.28 2.59 4.37 6.60 9.30 12.46 16.08 20.16 24.7
	n c (×10 5 ) 0.16 0.32 0.55 0.83 1.16 1.56	2.01	2.52 3.09
	nnz(A 0 )	5	5	5	5	5	5	5	5	5
	nnz(A 1 )	77	78	79	79	79	80	80	80	80
	nnz(A 2 )	299	315	316	322	332	335	337	339	341
	nnz(A 3 )	322	342	345	359	376	380	385	386	386
	nnz( P1 )	20	20	20	20	20	20	20	20	20
	nnz( P2 )	63	64	65	66	66	66	67	67	67
	nnz( P3 )	168	178	178	183	188	190	191	192	193

* This work was funded by CEA.


S Val Ère

Bitseki Penda

S Valère

KERNEL ESTIMATION OF THE TRANSITION DENSITY IN BIFURCATING MARKOV CHAINS

Keywords: Kernel estimator, cross validation method, rule of thumb type method, bifurcating Markov chains, binary trees, asymptotic normality Mathematics Subject Classification (2020): 62G05, 62G07, 62G20, 60J80, 60F05

We study the kernel estimator of the transition density of bifurcating Markov chains. Under some ergodic and regularity properties, we prove that this estimator is consistent and asymptotically normal. Next, in the numerical studies, we propose two data-driven methods to choose the bandwidth parameters. These methods are based on the so-called two bandwidths approach.

Introduction

This article is devoted to the study of the kernel estimators of the transition probability of bifurcating Markov chains. Before defining these estimators, let us first introduce useful definitions, notations and assumptions.

1.1. Bifurcating Markov chains. Let d ≥ 1 be a natural integer. In order to simplify the notations in the sequel, we set S = R d and we equip S with its Borel σ-algebra that we denote by S . We denote by B(S) (resp. B b (S), resp. B + (S)) the set of (resp. bounded, resp. non-negative) R-valued measurable functions defined on S. For f ∈ B(S), we set f ∞ = sup{|f (x)|, x ∈ S}. For a finite measure λ on (S, S ) and f ∈ B(S) we shall write λ, f for f (x) dλ(x) whenever this integral is well defined. We denote by C b (S) (resp. C + (S)) the set of bounded (resp. nonnegative) R-valued continuous functions defined on S. For all natural integer q ≥ 1, we equip S q with S ⊗q = S ⊗ . . . ⊗ S , the usual product σ-field on S q .

Let Q be a probability kernel on S × S , that is: Q(•, A) is measurable for all A ∈ S , and Q(x, •) is a probability measure on (S, S ) for all x ∈ S. For any f ∈ B b (S), we set for x ∈ S:

(1) (Qf )(x) = S f (y) Q(x, dy).
We define (Qf ), or simply Qf , for f ∈ B(S) as soon as the integral (1) is well defined, and we have Qf ∈ B(S). For n ∈ N, we denote by Q n the n-th iterate of Q defined by Q 0 = I d , the identity map on B(S), and

Q n+1 f = Q n (Qf ) for f ∈ B b (S).
Let P be a probability kernel on S × S ⊗2 , that is: P (•, A) is measurable for all A ∈ S ⊗2 , and P (x, •) is a probability measure on (S 2 , S ⊗2 ) for all x ∈ S. For any g ∈ B b (S 3 ) and h ∈ B b (S 2 ), we set for x ∈ S:

(2) (P g)(x) = S 2 g(x, y, z) P (x, dy, dz) and (P h)(x) = S 2 h(y, z) P (x, dy, dz).

1

We define (P g) (resp. (P h)), or simply P g for g ∈ B(S 3 )(resp. P h for h ∈ B(S 2 )), as soon as the corresponding integral (2) is well defined, and we have that P g and P h belong to B(S).

We now introduce some notations related to the regular binary tree. Recall that N is the set of non-negative integers and N * = N \ {0}. We set T 0 = G 0 = {∅}, G k = {0, 1} k and T k = 0≤r≤k G r for k ∈ N * , and T = r∈N G r . The set G k corresponds to the k-th generation, T k to the tree up the k-th generation, and T the complete binary tree. For i ∈ T, we denote by |i| the generation of i (|i| = k if and only if i ∈ G k ) and iA = {ij; j ∈ A} for A ⊂ T, where ij is the concatenation of the two sequences i, j ∈ T, with the convention that ∅i = i∅ = i.

We recall the definition of bifurcating Markov chain (BMC) from Guyon [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF].

Definition 1.1. We say a stochastic process indexed by T, X = (X i , i ∈ T), is a bifurcating Markov chain on a measurable space (S, S ) with initial probability distribution ν on (S, S ) and probability kernel P on S × S ⊗2 , a BMC in short, if: -(Initial distribution.) The random variable X ∅ is distributed as ν.

-(Branching Markov property.) For a sequence (g i , i ∈ T) of functions belonging to B b (S 3 ), we have for all k ≥ 0,

E i∈G k g i (X i , X i0 , X i1 )|σ(X j ; j ∈ T k ) = i∈G k Pg i (X i ).
We define three probability kernels P 0 , P 1 and Q on S × S by: P 0 (x, A) = P(x, A × S), P 1 (x, A) = P(x, S × A) for (x, A) ∈ S × S , and Q = 1 2 (P 0 + P 1 ).

Notice that P 0 (resp. P 1 ) is the restriction of the first (resp. second) marginal of P to S. Following Guyon [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF], we introduce an auxiliary Markov chain Y = (Y n , n ∈ N) on (S, S ) with Y 0 distributed as X ∅ and transition kernel Q. The distribution of Y n corresponds to the distribution of X I , where I is chosen independently from X and uniformly at random in generation G n . We shall write E x when X ∅ = x (i.e. the initial distribution ν is the Dirac mass at x ∈ S).

Let i, j ∈ T. We write i j if j ∈ iT. We denote by i ∧ j the most recent common ancestor of i and j, which is defined as the only u ∈ T such that if v ∈ T and v i, v j then v u. We also define the lexicographic order i ≤ j if either i j or v0 i and v1 j for v = i ∧ j. Let X = (X i , i ∈ T) be a BM C with kernel P and initial measure ν. For i ∈ T, we define the σ-field:

F i = σ(X u ; u ∈ T such that u ≤ i).
By construction, the σ-fields (F i ; i ∈ T) are nested as F i ⊂ F j for i ≤ j.

For i ∈ T and k ∈ N, we also define the σ-field:

H i,k = σ(X u , u ∈ iT k )
We end this section with a useful notations. By convention, for f, g ∈ B(S), we define the function f ⊗ g by (f ⊗ g)(x, y) = f (x)g(y) for x, y ∈ S and

f ⊗ sym g = 1 2 (f ⊗ g + g ⊗ f ) and f ⊗ 2 = f ⊗ f. Notice that P(g ⊗ sym 1) = Q(g) for g ∈ B + (S).
For all u ∈ T, we denote by X u = (X u , X u0 , X u1 ) the mother-daughters triangle. For a finite subset A ⊂ T, we define:

M A (f ) = u∈A f (X u ) if f ∈ B(S) and M A (f ) = u∈A f (X u ) if f ∈ B(S 3 ).
In the sequel we will also use the following notation: let g and h be two functions which depend on one variable, x say; we denote by g ⊕ h the function of three variables, xx 0 x 1 := (x, x 0 , x 1 ) say, defined by (g ⊕ h)(xx 0 x 1 ) = g(x 0 ) + h(x 1 ).

1.2. Assumptions on the law of the bifurcating Markov chains (X i , i ∈ T). For a set F ⊂ B(S) of R-valued functions, we write

F 2 = {f 2 ; f ∈ F }, F ⊗ F = {f 0 ⊗ f 1 ; f 0 , f 1 ∈ F },
and P (F ) = {P f ; f ∈ F } whenever a kernel P act on F . Following [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF], we state a structural assumption on the set of functions we shall consider.

Assumption 1.2. Let F ⊂ B(S) be a set of R-valued functions such that: (i) F is a vector subspace which contains the constants;

(ii) F 2 ⊂ F ; (iii) F ⊂ L 1 (ν); (iv) F ⊗ F ⊂ L 1 (P(x, •
)) for all x ∈ S, and P(F ⊗ F ) ⊂ F .

The condition (iv) implies that P 0 (F ) ⊂ F ,P 1 (F ) ⊂ F as well as Q(F ) ⊂ F . Notice that if f ∈ F , then even if |f | does not belong to F , using conditions (i) and (ii), we get, with g = (1+f 2 )/2, that |f | ≤ g and g ∈ F . Typically, the set F can be the set C b (S) of bounded real-valued functions, or the set of smooth real-valued functions such that all derivatives have at most polynomials growth.

Following [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF], we also consider the following ergodic properties for Q. Assumption 1.3. There exists a probability measure µ on (S, S ) such that F ⊂ L 1 (µ) and for all f ∈ F , we have the point-wise convergence lim n→∞ Q n f = µ, f and there exists g ∈ F with:

(3)

|Q n (f )| ≤ g for all n ∈ N.

Moreover, there exists a function V : [1, +∞) → (0, ∞) such that V ∈ F and constants α ∈ (0, 1) and M < ∞ such that:

(4) sup |f |≤V |Q n f -µ, f | ≤ M α n V for all n ∈ N.
Remark 1.4. In particular, (4) implies that for all f ∈ B b (S), we have

(5) |Q n f -µ, f | ≤ M f ∞ α n V for all n ∈ N.
Next, we have the following assumption on the existence of the density of P.

Assumption 1.5. The transition kernel P has a density, still denoted by P, with respect to the Lebesgue measure.

Remark 1.6. Assumption 1.5 implies that the transition kernel Q has a density, still denoted by Q, with respect to the Lebesgue measure. More precisely, we have Q(x, y) = 2 -1 S (P(x, y, z) + P(x, z, y))dz. This implies in particular that the invariant probability µ has a density, still denoted by µ, with respect to the Lebesgue measure (for more details, we refer for e.g. to [START_REF] Duflo | Random iterative models[END_REF], chap 6).

Remark 1.7. Under Assumption 1.5, the probability measure µ defined on S 3 by µ (dxx 0 x 1 ) = µ(dx)P (x, dx 0 , dx 1 ), has density with respect to the Lebesgue measure, that we also denote by µ , given by µ (xx 0 x 1 ) = µ(x)P(x, x 0 , x 1 ), for all xx 0 x 1 ∈ S 3 . Assumption 1.8. We assume that the following constant is finite:

C 0 = sup x,x0,x1∈S (µ(x) + Q(x, x 0 ) + P(x, x 0 , x 1 )).
Remark 1.9. We recall from [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF]Theorem 11 and Corollary 15] that under Assumptions 1.2 and 1.3, we have for f ∈ F the following convergence in L 2 (resp. a.s.): [START_REF] Comte | On a nadaraya-watson estimator with two bandwidths[END_REF] lim

n→∞ |G n | -1 M Gn (f ) = µ, f and lim n→∞ |T n | -1 M Tn (f ) = µ, f .
Now, the rest of the paper is organized follows. In Section 2, we define the estimators of the transition density P based on the observation of a subpopulation. We will see that these are quotient estimators. In Section 3, we study the consistency and the asymptotic normality of the numerators of the estimators of P. Section 4 is dedicated to the study of consistency and asymptotic normality of the estimators of P. In Section 5, we will illustrate the consistency of our estimators in a bifurcating Markov model called bifurcating autoregressive process (BAR, for short). In particular, we will develop two data-driven bandwidth selection methods: the least squares Cross-Validation in Section 5.1 and the rule of thumb type method in Section 5.2. Sections 6-8 are dedicated to the proofs of the main Theorems. In Section 9, we prove a useful inequality and in Section 10, we recall some useful results.

Kernel estimators of the transition density P

Recall that S = R d . Our aim is to estimate the transition density P from the observation of the subpopulation A n ∈ {G n , T n }. For that purpose, assume we observe X n = (X u ) u∈An i.e. we have 2 n+2 -1 (or 3 × 2 n ) random variables with value in S. Let K 0 : S → R and K : S 3 → R be a functions such that S K 0 (x)dx = 1 and K = K 0 ⊗K 0 ⊗K 0 . We also have S 3 K(xx 0 x 1 )dxx 0 x 1 = 1. Let (h n , n ∈ N) be a sequence of positive numbers which converges to 0 as n goes to infinity. When there is no ambiguity, we write h for h n . Let A n ∈ {T n , G n }. We define, for all x ∈ S:

(7) µ An (x) = 1 |A n |h d/2 u∈An K 0hn (x -X u ), where K 0hn (x -y) = h -d/2 n K 0 (h -1 n (x -y))
and for all xx 0 x 1 ∈ S 3 :

(8) µ An (xx 0 x 1 ) = 1 |A n |h 3d/2 u∈An K hn (xx 0 x 1 -X u ) and P An (xx 0 x 1 ) = µ An (xx 0 x 1 ) µ An (x) , where K hn (xx 0 x 1 -yy 0 y 1 ) = h -3d/2 n K h -1 n (x -y), h -1 n (x 0 -y 0 ), h -1 n (x 1 -y 1 )
, with the convention that P An (xx 0 x 1 ) = 0 if µ An (x) = 0. However, we stress that if we assume that K 0 is strictly positive, then µ An (x) > 0 for all x ∈ S.

From now on, we fix xx 0 x 1 ∈ S 3 , that is, we are interested in the estimation at the point xx 0 x 1 . We assume that µ(x) = 0. We consider the function f n defined by: [START_REF] Gyori | Hypothesis testing for markov chain monte carlo[END_REF] f n (yy 0 y 1 ) = K h (xx 0 x 1yy 0 y 1 ).

If we want to be more rigorous, we must write f n,xx0x1 instead of f n . But, we choose to write without the index xx 0 x 1 in order to simplify the writing.

Remark 2.1. Note that asymptotic behavior (consistence and asymptotic normality) of µ An have been studied in [START_REF] Bitseki Penda | Central limit theorem for bifurcating markov chains[END_REF].

Remark 2.2. We stress that the results of this paper can be straightforward extended the case where the bandwidth h h h is a vector of R 3d , with possibly different coordinates. More precisely, one can take the bandwidth h h h = (h i , 1 ≤ i ≤ 3d), where the h i 's may take different values. For our convenience, we choose to work with the case where all the coordinates are the same, that is h i = h for all 1 ≤ i ≤ 3d.

3. Consistency and Asymptotic normality for µ An (xx 0 x 1 )

First, we will study the consistency and the asymptotic normality of µ An (xx 0 x 1 ). We set fn = f nµ, Pf n . We begin with the study asymptotic normality of

N n,∅ (f n ) = |G n | -1/2 M An ( fn ).
This is motivated by the following decomposition: [START_REF] Hall | Martingale limit theory and its application[END_REF]

µ An (xx 0 x 1 )-µ (xx 0 x 1 ) = (|A n | |G n | -1/2 h 3d/2 ) -1 N n,∅ (f n ) + (h -3d/2 µ , f n -µ (xx 0 x 1 )).
We will need the following assumption on the bandwidth and on the kernel. Assumption 3.1.

We assume that: (i) h n = 2 -nγ and 2α 2 < 2 3dγ for some γ ∈ (0, 1/3d).

(ii) The kernel K 0 (resp. K 2 0 ) is integrable and square integrable. Remark 3.2. Assumption 3.1, (i) implies in particular that [START_REF] Parzen | On estimation of a probability density function and mode[END_REF] lim

n→∞ |G n |h 3d n = ∞ and lim n→∞ (2α 2 ) n h 3d n = 0.
Note that Assumption 3.1, (i) is automatically satisfied if 2α 2 ≤ 1, regardless of the value of α. For 2α 2 > 1, this Assumption implies that the choice of the bandwidth is function of the ergodicity rate of the auxiliary Markov chain Y .

We have the following result.

Theorem 3.3. Let X be a BMC with kernel P and initial distribution ν such that Assumptions 1.2, 1.3, 1.5, 1.8 and 3.1 hold. Then, we have the following convergence in distribution:

N n,∅ (f n ) (d) ----→ n→∞ G,
where G is a centered Gaussian random variable with finite variance

σ 2 = 2 K 0 6 2 µ (x, x 0 , x 1 ) if A n = T n and σ 2 = K 0 6 2 µ (x, x 0 , x 1 ) if A n = G n .
Proof. The proof of Theorem 3.3 is postponed to Section 6.

Next, in order to study the asymptotic normality of µ An (xx 0 x 1 ), we do the following additional hypothesis.

Assumption 3.4. We assume that Assumption 3.1 holds and there exists s > 0 such that the following holds.

(iv) The density µ (resp. µ) belongs to the (isotropic) Hölder class of order (s, . . . , s) ∈ R 3d (resp. (s, . . . , s) ∈ R d ): The density µ admits partial derivatives with respect to x j , for all j ∈ {1, . . . 3d}, up to the order s and there exists a finite constant L > 0 such that for all x = (x 1 , . . . , x 3d ), ∈ R 3d , t ∈ R and j ∈ {1, . . . , 3d}:

∂ s µ ∂x s j (x -j , t) - ∂ s µ ∂x s j (x) ≤ L|x j -t| {s} ,
where (x -j , t) denotes the vector x where we have replaced the j th coordinate x j by t, with the convention ∂ 0 µ /∂x 0 j = µ . The same thing for the density µ.

(v) The kernel K 0 is of order ( s , . . . , s ) ∈ N d : We have R d |x| s K 0 (x) dx < ∞ and R x k j K 0 (x) dx j = 0 for all k ∈ {1, . . . , s } and j ∈ {1, . . . , d}. (vi) Bandwith control: We have γ > 1/(2s + 3d), that is lim n→∞ |G n |h 2s+3d n = 0.
Notice that Assumption 3.4-(iv) implies that µ (resp. µ) is at least Hölder continuous as s > 0. We have the following result.

Theorem 3.5. Let X be a BMC with kernel P and initial distribution ν. Under Assumptions of Theorem 3.3, we have for all (x, x 0 , x 1 ) and

A n ∈ {G n , T n } (12) µ An (xx 0 x 1 ) P ----→ n→∞ µ (xx 0 x 1 ) in probability.
Moreover, under the additional Assumption 3.4, we have the following convergence in distribution:

|A n | 1/2 h 3d/2 n ( µ An (xx 0 x 1 ) -µ (xx 0 x 1 )) (d) ----→ n→∞ G in distribution,
where G is a centered Gaussian random variable with finite variance σ 2 = K 0 6 2 µ (x, x 0 , x 1 ). Proof. The proof is postponed to Section 7.

Consistency and Asymptotic normality for

P An (xx 0 x 1 )
We are now in position to state consistency and asymptotic normality of kernel estimator of the transition density P. First, as a consequence of ( 8), ( 12) and (43) below, we have the following result.

Lemma 4.1. Under the Assumptions of Theorem 3.3, we have for all (x, x 0 , x 1 ) and A n ∈ {G n , T n } :

P An (x, x 0 , x 1 ) P ----→ n→∞ P(x, x 0 , x 1 ) in probability.
Next, we have the following result. Theorem 4.2. Let X be a BMC with kernel P and initial distribution ν. Under the assumptions of Theorem 3.3 and the additional Assumption 3.4, we have,

|A n |h 3d n ( P An (x, x 0 , x 1 ) -P(x, x 0 , x 1 )) (d) ----→ n→∞ G,
where G is a centered Gaussian real-valued random variable with mean 0 and variance

σ 2 = K 0 6 2 P(x, x 0 , x 1 )/µ(x).
Proof. The proof is postponed to Section 8.

Numerical studies

We consider the real-valued Gaussian bifurcating autoregressive process (BAR) X = (X u , u ∈ T) where X ∅ is arbitrary and for all u ∈ T:

(13) X u0 = a 0 X u + b 0 + ε u0 X u1 = a 1 X u + b 1 + ε u1 , with a 0 , a 1 ∈ [-1, 1], b 0 , b 1 ∈ R and ((ε u0 , ε u1 ), u ∈ T)
an independent sequence of bivariate Gaussian N(0, Γ) random vectors independent of X ∅ with covariance matrix, with σ > 0 and ρ ∈ R such that |ρ| < σ 2 :

Γ = σ 2 ρ ρ σ 2 .
Then the process X = (X u , u ∈ T) is a BMC with transition probability P given by:

P(x, dy, dz) = 1 2π σ 4 -ρ 2 exp - σ 2 2(σ 4 -ρ 2 ) g(x, y, z) dydz, with g(x, y, z) = (y -a 0 x -b 0 ) 2 -2ρσ -2 (y -a 0 x -b 0 )(z -a 1 x -b 1 ) + (z -a 1 x -b 1 ) 2 .
The transition kernel Q of the auxiliary Markov chain is defined by:

Q(x, dy) = 1 2 √ 2πσ 2 e -(y-a0x-b0) 2 /2σ 2 + e -(y-a1x-b1) 2 /2σ 2 dy.
We will estimate the transition density P in a compact set D ⊂ R 3 . For that purpose, we use the estimator P Gn (xx 0 x 1 ), for all xx 0 x 1 ∈ D, given in ( 8), with the Gaussian kernel K 0 defined by

(14) K 0 (x) = 1 √ 2π e -x 2 /2 .
Since the bandwidth is a function of the ergodicity rate which is unknown, we have to develop a method based on data in order to select it. To select the optimal bandwidth for P Gn defined in ( 8), we will use the so-called "two bandwidths approach" (see for e.g. [START_REF] Comte | On a nadaraya-watson estimator with two bandwidths[END_REF]). More precisely, since P Gn is a quotient estimator, we select separately the bandwidths for the numerator (h N , say) and the denominator (h D , say). For that purpose, we propose two methods: the cross validation and the rule of thumb type method. The objective here is not to study nor to compare theoretically these two methods. This will be done in the future works. Our objective is only the see the numerical performances of each method. Our conclusion is that even if the rule of thumb developed in this paper give a crude approximation, it as more computational benefit with respect to the least squared cross validation.

5.1.

Bandwidth selection by least squares Cross-Validation method. We choose the bandwidths which minimises the mean integrated squared errors (MISEs)

E R 3 ( µ Gn -µ ) 2 (xx 0 x 1 )dxdx 0 dx 1 and E R ( µ Gn -µ) 2 (x)dx ,
where µ Gn and µ Gn are defined in ( 7) and ( 8). This is equivalent to minimise the functions J and J defined by

J (h) = E R 3 ( µ Gn ) 2 (xx 0 x 1 )dxdx 0 dx 1 -2E R 3 ( µ Gn µ )(x, x 0 , x 1 )dxdx 0 dx 1 and J(h) = E R ( µ Gn ) 2 (x)dx -2E R ( µ Gn µ)(x)dx .
The method to select the bandwidths is the following.

(1) We divide the sample (X u )

u∈Gn into K disjoints subsamples {(X u ) u∈G (k) n , k ∈ {1, . . . , K}}, with (G (k) n , k ∈ {1, . . . , K}) a partition of G n . (2) For each subsample (X u ) u∈G (k) n : (a)
We set µ Gn[-k] and µ Gn[-k] the estimators of µ and µ obtaining using the subsample

(X u ) u∈Gn K (X u ) u∈G (k) n
, where for two sets A and B, B K A denotes the set of elements in B but not in A. More precisely,

µ Gn[-k] (xx 0 x 1 ) = 1 |G n [-k]|h 3 u∈Gn[-k] K 0 x -X u h K 0 x 0 -X u0 h K 0 x 1 -X u1 h and µ Gn[-k] (x) = 1 |G n [-k]|h u∈Gn[-k] K 0 x -X u h ,
where we set

G n [-k] = G n K G (k)
n . (b) We approximate J and J by

J (k) (h) = R ( µ Gn[-k] ) 2 (x)dx - 2 |G (k) n | u∈G (k) n µ Gn[-k] (X u ) and
J (k) (h) = R 3 ( µ Gn[-k] ) 2 (xx 0 x 1 )dxdx 0 dx 1 - 2 |G (k) n | u∈G (k) n µ Gn[-k] (X u , X u0 , X u1 ).
(3) Let H = {h 1 , . . . , h m } ⊂ (0, 1] be a bandwidth grid. Then, the selected bandwidths h N and h D for the numerator and the denominator of P Gn are given by:

h N := arg min h∈H 1 K K k=1 J (k) (h) and h D := arg min h∈H 1 K K k=1 J (k) (h).
Finally, the estimator used for numerical studies is P Gn defined by:

P Gn (xx 0 x 1 ) = µ An (xx 0 x 1 ) µ An (x) , with µ Gn (x) = 1 |G n | h D u∈Gn K 0 x -X u h D and µ Gn (xx 0 x 1 ) = 1 |G n | h 3 N u∈Gn K 0 x -X u h N K 0 x 0 -X u0 h N K 0 x 1 -X u1 h N .
This method is known as the K-fold cross validation. One advantage of this method in the context of bifurcating Markov chains is that it is not requires the knowledge of the ergodicity rate. The main drawback being that it requires a lot of time for calculations.

Gaussian symmetric BAR reference bandwidth selection.

In order to define a selection rule, we consider the special case of Gaussian BAR defined by (13) where a 0 = a 1 := a and ρ = 0 as a reference model. It is well known (see [START_REF] Bitseki Penda | Central limit theorem for bifurcating markov chains under pointwise ergodic conditions[END_REF]) that the densities of the transition kernel Q of the auxiliary Markov chain and the invariant probability µ associated to Q are given by:

Q(x, y) = 1 √ 2πσ 2 exp - (y -ax) 2 2σ 2 and µ(x) = 1 √ 2πσ a exp - x 2 2σ 2 a ,
where

σ a = σ/ √ 1 -a 2 .
The density of the transition probability P associated to this bifurcating Markov chain is defined by P(x, y, z) = Q(x, y)Q(x, z) and then, µ (x, y, z) = µ(x)Q(x, y)Q(x, z). We then have that the invariant densities µ and µ are square integrable and twice differentiable, with the second order derivative of µ and all the second order partial derivatives of µ bounded, continuous and square integrable. It is also well that the Markov chain with transition Q is geometrically ergodic and that the geometric ergodic rate of convergence is a (for more details, see for e.g Example 2.8 in [START_REF] Bitseki Penda | Central limit theorem for bifurcating markov chains under pointwise ergodic conditions[END_REF]). In particular, following the proof of Proposition 28 in [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF], one can prove that for all derivable function f such f and f are bounded, we have

(15) |Q n f (x) -µ, f | ≤ f ∞ (σ(1 -a) -1 + |x|)a n .
We assume that L(X ∅ ) = µ, that is X ∅ is distributed as µ, which implies that the process is stationary. We are now going to behave as if we did not know the invariant measures µ and µ and the transition probability P. Recall the kernel density estimator of µ defined in [START_REF] Duflo | Random iterative models[END_REF] and the kernel K 0 defined in (14). Recall also the kernel estimator of the transition density P defined in [START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF]. To ease notation, we write µ and µ instead of µ Gn and µ Gn respectively. We recall that our strategy is to select bandwidth for the numerator and the denominator in the estimation of P.

First, we treat the denominator µ. The selection rule is based on the following asymptotic upper bound, known as asymptotic mean squared error:

(16) E ( µ(x) -µ(x)) 2 ≤ h 4 2 κ 2 2 µ (x) 2 + 4 K 0 2 2 |G n |h µ(x) + 2 C a,σ a 2 |G n | n-1 k=1 (2a 2 ) k + O(1).
where

(17) C a,σ = 1 πσ 2 (1 + a)(1 -a) 2 and κ 2 = R y 2 K 0 (y)dy = 1.
We postponed the proof of (16) in Section 9. Now, let p be a non negative probability density defined in R such that p ∞ ≤ 1. Then, (16) implies that

(18) R E ( µ(x) -µ(x)) 2 p(x)dx ≤ h 4 2 κ 2 2 R µ (x) 2 dx + 4 K 0 2 2 |G n |h + 2 C a,σ a 2 |G n | n-1 k=1 (2a 2 ) k + O(1).
The term in the left hand side of (18) is a modification of asymptotic mean integrated squared error that we call p-AMISE. We have introduced it because the last term in (16) does not depend on x. Finally, (18) suggests us to choose the bandwidth which minimises the function G defined by

G(h) = h 4 2 κ 2 2 R µ (x) 2 dx + 4 K 0 2 2 |G n |h + M a,σ a 2n 1 {2a 2 >1} ,
where M a,σ = 2/(πσ 2 a 2 (1+a)(1-a) 2 (2a 2 -1)). Optimizing in h, we get that the optimal bandwidth h D (for the denominator of P Gn defined in ( 8)) is given by

h D = (c 1 /4c 2 ) 1/5 |G n | -1/5 1 {2a 2 ≤2 -1/5 } + (c 1 /M a,σ ) (2a 2 ) -n 1 {2a 2 >2 -1/5 } ,
where

c 1 = 4 K 0 2 2 = 2 √ π and c 2 = (1/2)κ 2 2 R µ (x) 2 dx = 3 16 √ π σ -5 a .
Next, we treat the numerator µ of P Gn . Recalling Remark 2.2, we consider the general case where for all xx 0 x 1 ∈ S 3 :

µ Gn (xx 0 x 1 ) = 1 |G n |hh 0 h 1 u∈Gn K 0 (h -1 (x -X u ))K 0 (h -1 0 (x 0 -X u0 ))K 0 (h -1 1 (x 1 -X u1 )).
Recall that for a vector v v v, v v v t denotes its transpose. As in (16), we have the following asymptotic upper bound:

(19) E µ (xx 0 x 1 ) -µ (xx 0 x 1 ) 2 ≤ 1 2 κ 2 2 h h h t H H H(xx 0 x 1 ) h h h 2 + 6 K 0 6 2 |G n |hh 0 h 1 µ (xx 0 x 1 ) + C a,σ P(xx 0 x 1 ) 2 |G n | n-1 k=1 (2 a 2 ) k + O(1),
where κ 2 is defined in (17), h h h = (h, h 0 , h 1 ) t , and

H H H(xx 0 x 1 ) =     ∂ 2 µ ∂x 2 (x, x 0 , x 1 ) 0 0 0 ∂ 2 µ ∂x 2 0 (x, x 0 , x 1 ) 0 0 0 ∂ 2 µ ∂x 2 1 (x, x 0 , x 1 )     , C a,σ = 4 e πσ 2 a 2 (1 -a) 2 (1 + a)
.

We let the proof of (19) to the reader since it follows the same lines that of (16). Let p be a non negative probability density defined in R such that p ∞ ≤ 1. Integrating (19) with respect to p(x)dxdx 0 dx 1 , we get

R 3 E µ (xx 0 x 1 ) -µ (xx 0 x 1 ) 2 p(x)dxdx 0 dx 1 ≤ 1 2 κ 2 2 R 3 h h h t H H H(xx 0 x 1 ) h h h 2 dxdx 0 dx 1 + 6 K 0 6 2 |G n |hh 0 h 1 + C a,σ 4πσ 2 |G n | n-1 k=1 (2 a 2 ) k + O(1).

Now, the latter equation suggests us to choose the vector bandwidth h h h which minimises the function G defined by

G (h, h 0 , h 1 ) = 1 2 κ 2 2 R 3 h h h t H H H(xx 0 x 1 ) h h h 2 dxdx 0 dx 1 + 6 K 0 6 2 |G n |hh 0 h 1 + M a,σ a 2n 1 {2a 2 >1} .
where M a,σ = (1/(4 π σ 2 (2a 2 -1))) C a,σ . Optimizing the function G in h h h, we get that the optimal bandwidth h h h N = (h N , h 0N , h 1N ) is given by

h N = c 2 /(4c ) 1/7 |G n | -1/7 1 {2a 2 ≤2 -3/7 } + (c 2 /M a,σ ) 1/3 (2a 2 ) -n/3 1 {2a 2 >2 -3/7 } , h 0N = c 2 /(4c 0 ) 1/7 |G n | -1/7 1 {2a 2 ≤2 -3/7 } + (c 2 /M a,σ ) 1/3 (2a 2 ) -n/3 1 {2a 2 >2 -3/7 } , h 1N = c 2 /(4c 1 ) 1/7 |G n | -1/7 1 {2a 2 ≤2 -3/7 } + (c 2 /M a,σ ) 1/3 (2a 2 ) -n/3 1 {2a 2 >2 -3/7 } , where κ 2 is defined in (17), c 2 = 6 K 0 6 2 = 6/(8π √ π) and c = 1 2 κ 2 2 R 3 ∂ 2 µ ∂x 2 (x, x 0 , x 1 ) 2 dxdx 0 dx 1 = 3(1 + a 2 ) 2 64π √ π(1 -a 2 ) 3 σ -7 a , c 0 = 1 2 κ 2 2 R 3 ∂ 2 µ ∂x 2 0 (x, x 0 , x 1 ) 2 dxdx 0 dx 1 = 3 64π √ π(1 -a 2 ) 3 σ -7 a , c 1 = 1 2 κ 2 2 R 3 ∂ 2 µ ∂x 2 1 (x, x 0 , x 1 ) 2 dxdx 0 dx 1 = 3 64π √ π(1 -a 2 ) 3 σ -7 a .
We have

c 1 4c 2 = b 1 σ 5 a , c 1 M a,σ = b 2 σ 2 a , c 2 4 c = b 3 σ 7 a , c 2 4 c 0 = c 2 4 c 1 = b 4 σ 7 a , c 2 M a,σ = b 5 σ 4 a ,
where

b 1 = 32 3 , b 2 = √ πa 2 (1 -a 2 )(1 + a)(1 -a) 2 (2a 2 -1), b 3 = 48(1 -a 2 ) 3 12(1 + a 2 ) 2 , b 4 = 4(1 -a 2 ) 3 , b 5 = 3a 2 (1 -a 2 ) 3 (1 + a)(2a 2 -1) 4 √ π .
Since for a ∈ (0, 1) the constants b i , i ∈ {1, . . . , 5}, are bounded, we can approximate h D , h N , h 0N , h 1N by:

h D = |G n | -1/5 σ a 1 {2 a 2 ≤2 -1/5 } + (2 a 2 ) -n σ a 1 {2 a 2 >2 -1/5 } , h N = h 0N = h 1N = |G n | -1/7 σ a 1 {2 a 2 ≤2 -3/7 } + (2 a 2 ) -n/3 σ a 1 {2 a 2 >2 -3/7 } ,
where σ a is the estimator of the standard deviation of the measure µ and a is the estimator of the geometric ergodic rate. Note that in practice, the estimators h N , h 0N and h 1N differ slightly. Indeed, for h D and h N , σ a is computed using the sample (X u , u ∈ G n ), for h 0N , σ a is computed using the sample (X u0 , u ∈ G n ) and for h 1N , σ a is computed using the sample (X u1 , u ∈ G n ).

Recall that for i ∈ T and A ⊂ T, iA = {ij, j ∈ A}, where ij is the concatenation of the two sequences i, j ∈ T. For the geometric ergodic rate, we propose the following estimator, which is inspired from [START_REF] Gyori | Hypothesis testing for markov chain monte carlo[END_REF]:

a = u∈Gn-m+1 v∈uGm-1 (X u -X)(X v -X) u∈Gn (X u -X) 2 1/m with X = 1 |G n | u∈Gn X u ,
where m is a large enough natural integer such that m = O(n). The choice m = n/2 + 1 seems to be relevant. Finally, the estimator used for numerical studies is P Gn defined by:

P Gn (xx 0 x 1 ) = µ An (xx 0 x 1 ) µ An (x) , with µ Gn (x) = 1 |G n | h D u∈Gn K 0 x -X u h D and µ Gn (xx 0 x 1 ) = 1 |G n | h 3 N u∈Gn K 0 x -X u h N K 0 x 0 -X u0 h 0N K 0 x 1 -X u1 h 1N .
This method is an adaptation of the rule of thumb developed by Silverman in [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF]. The novelty here is that the ergodic rate of convergence is taken into account in the estimation procedure.

In the context of BMC, The main advantage of this method is that it not requires a lot of time for calculations. However, this method is a crude approximation which works for approximately "Gaussian" bifurcating Markov chains.

Numerical illustrations.

In order to validate our method, we consider two cases: case 1: (a 0 , a 1 , b 0 , 1 , σ, ρ) = (0.7, 0.5, 0, 0, 1, 0); case 2: (a 0 , a 1 , b 0 , 1 , σ, ρ) = (1.2, 0.7, 0, 0, 1, 0); In case 2, we allow the dynamic of the new pole to be unstable, even if the entire dynamic of the system is stable. Following the terminology of Bitseki and Delmas in [START_REF] Bitseki Penda | Central limit theorem for bifurcating Markov chains under L 2 ergodic conditions[END_REF][START_REF] Bitseki Penda | Central limit theorem for bifurcating markov chains under pointwise ergodic conditions[END_REF], the case 2 corresponds to supercritical case.

As we can see, Figure 1-8, the two methods allow to recover the true function when the size of the data increases. Consequently, we conclude that our method is valid. We begin the proof with A n = T n . Let (p n , n ∈ N) be a non-decreasing sequence of elements of N * such that, for all λ > 0:

p n < n, lim n→∞ p n /n = 1 and lim n→∞ n -p n -λ log(n) = +∞.
When there is no ambiguity, we write p for p n . Recall the function f n defined in [START_REF] Gyori | Hypothesis testing for markov chain monte carlo[END_REF]. We have the following decomposition:

(20) N n,∅ (f n ) = R 0 (n) + R 1 (n) + ∆ n (f n ),
where: The projection of PGn (black) and P (green) on the line x = 3; x1 = 0. We use the rule of thumb type method and n = 15; a0 = 1.2; a1 = 0.7; x0 P (black) and P (green)

R 0 (n) = |G n | -1/2 n-p-1 k=0 M G k ( fn ); R 1 = i∈Gn-p E[N n,i (f n )|F i ]; ∆ n (f n ) = i∈Gn-p ∆ n,i (f n ),
and for all i ∈ G n-p ,

N n,i (f n ) = |G n | -1/2 p =0 M iG p-( fn ) and ∆ n,i (f n ) = N n,i (f n ) -E[N n,i (f n )|F i ].
Note that using the branching Markov property, we have, for all i ∈ G n-p , ( 21)

E[M iG p-( fn )|F i ] = E Xi [M G p-(P fn )].
We have the following convergence.

Lemma 6.1. Under the assumptions of Theorem 3.3, we have that

lim n→∞ E[R 0 (n) 2 ] = 0.
Proof. We have

E[R 0 (n) 2 ] = |G n | -1 E[( n-p-1 k=0 u∈G k fn (X u ) 2 ] ≤ |G n | -1 ( n-p-1 k=0 E[( u∈G k fn (X u ) 2 ] 1/2 ) 2 , ( 22
)
where we used the Minkowski inequality for the first inequality. By developing the term in the expectation, we get

E[( u∈G k fn (X u ) 2 ] = E[ u =v∈G k E[ fn (X u ) fn (X v )|X u , X v ]] + E[ u∈G k E[( fn ) 2 (X u )|X u ]] = E[ u =v∈G k P fn (X u )P fn (X v )] + E[ u∈G k P(( fn ) 2 )(X u )] = E[( u∈G k P fn (X u )) 2 ] + E[ u∈G k (P((f n ) 2 ) -(Pf n ) 2 )(X u )],
where we used the branching Markov property for the second inequality and the fact that P(( fn ) 2 )-(P fn ) 2 = P((f n ) 2 ) -(Pf n ) 2 for the third equality. Using (22) and using the inequalities

√ a + b ≤ √ a + √ b and (a + b) 2 ≤ 2a 2 + 2b 2 , we get E[R 0 (n) 2 ] ≤ |G n | -1 ( n-p-1 k=0 (E[( u∈G k P fn (X u )) 2 ] 1/2 + E[ u∈G k (P((f n ) 2 ) -(Pf n ) 2 )(X u )] 1/2 )) 2 ≤ 2|G n | -1 (( n-p-1 k=0 E[M G k (P( fn )) 2 ] 1/2 ) 2 + ( n-p-1 k=0 E[M G k (P((f n ) 2 ) -(Pf n ) 2 )] 1/2 ) 2 ).
Note that from Lemma 10.1, we have

P(f n ) ∞ ≤ Ch d/2 n and P(f 2 n ) ∞ ≤ C. Recall P( fn ) = P(f n ) -µ, P(f n ) .
Then, using (53), ( 5) and Lemma 10.1, we get

E[M G k (P fn ) 2 ] ≤ C h d n if k ∈ {0, 1} and for all k ≥ 2, E[M G k (P fn ) 2 ] ≤ C h 2d n 2 k if 2α 2 ≤ 1 C 2 k (h 2d n + (2α 2 ) k h 3d n ) if 2α 2 > 1.
It follows for the two last inequalities that

|G n | -1 ( n-p-1 k=0 E[M G k (P( fn )) 2 ] 1/2 ) 2 ≤ C 2 -n h 3d n + C 2 -p h 2d n + 2 -p (2α 2 ) n-p h 3d n 1 {2α 2 >1} .
Using [START_REF] Parzen | On estimation of a probability density function and mode[END_REF], Assumption 3.1 and since lim n→∞ p n = ∞, it follows that

lim n→∞ |G n | -1 ( n-p-1 k=0 E[M G k (P( fn )) 2 ] 1/2 ) 2 = 0.
Next, using Lemma 10.1, we get

E[M G k (P((f n ) 2 ) -(Pf n ) 2 )] ≤ C 2 k . This implies that lim n→∞ |G n | -1 ( n-p-1 k=0 E[M G k (P((f n ) 2 ) -(Pf n ) 2 )] 1/2 ) 2 ≤ C lim n→∞ 2 -p = 0
and this ends the proof.

Next, we have the following convergence.

Lemma 6.2. Under the assumptions of Theorem 3.3, we have that

lim n→∞ E[R 1 (n) 2 ] = 0.
Proof. Using (21), we get

R 1 (n) = p k=0 R 1 (k, n), with R 1 (k, n) = |G n | -1/2 |G p-k | M Gn-p (Q p-k P fn ).
It follows that

(23) E R 1 (n) 2 1/2 ≤ p k=0 E R 1 (k, n) 2 1/2 .
Following the proof of Lemma 4.2 in [3] and using (5) and Lemma 10.1, we find that

E R 1 (k, n) 2 ≤ C 2 -p h 2d n 1 {k=p} +      C h 3d n 2 -p (2α) 2(p-k) if 2α 2 < 1 C (n -p) h 3d n 2 -k if 2α 2 = 1 C 2 -p (2α 2 ) n-p h 3d n (2α) 2(p-k) if 2α 2 > 1. From (23), this implies that E R 1 (n) 2 1/2 ≤ C 2 -p/2 h d n +      C h 3d/2 n p k=0 2 -k/2 (2α 2 ) (p-k)/2 if 2α 2 < 1 C (n -p) 1/2 h 3d/2 n if 2α 2 = 1 C (2α 2 ) n/2 h 3d/2 n if 2α 2 > 1.
From the latter inequality and using Assumption 3.1, we deduce that lim

n→∞ E[R 1 (n) 2 ] = 0.
We now study the bracket

V (n) = i∈Gn-p E[∆ n,i (f n ) 2 |F i ].
Note that for i ∈ G n-p , we have

E[∆ n,i (f n ) 2 |F i ] = |G n | -1 E[( p k=0 M iG p-k ( fn )) 2 |F i ] -|G n | -1 (E[ p k=0 M iG p-k ( fn )|F i ]) 2 .
Using the branching Markov chain property, this implies that (24)

V (n) = V 1 (n) + V 2 (n) -R 2 (n), with V 1 (n) = |G n | -1 i∈Gn-p p k=0 E Xi [M G p-k ( fn ) 2 ], V 2 (n) = 2|G n | -1 i∈Gn-p 0≤k< ≤p E Xi [M G p-( fn )M G p-k ( fn )], (25) R 2 (n) = |G n | -1 i∈Gn-p ( p k=0 2 p-k Q p-k P fn (X u )) 2 .
We have the following result. Lemma 6.3. Under the assumptions of Theorem 3.3, we have the following convergence:

lim n→∞ E[R 2 (n)] = 0.
Proof. We have using (52), ( 5) and Lemma 10.1:

E [R 2 (n)] = |G n | -1 |G n-p | ν, Q n-p p k=0 |G p-k | Q p-k P fn 2 ≤ C 2 -p ν, Q n-p ((P fn ) 2 ) + C 2 -p ν, Q n-p (( p-1 k=0 2 p-k Q p-k-1 (QP fn )) 2 ) ≤ C 2 -p h 2d n + C 2 -p h 3d n a n
, where the sequence (a n , n ≥ 1) is defined by

a n =      1 if 2α < 1 p 2 if 2α = 1 (2α) 2p if 2α > 1.
Using Assumption 3.1, and in particular Remark 3.2, it follows that lim n→∞ E[R 2 (n)] = 0.

Next, we have the following result. Lemma 6.4. Under the assumptions of Theorem 3.3, we have the following convergence:

lim n→∞ E[V 2 (n) 2 ] = 0.
Proof. Let 0 ≤ k < ≤ p and i ∈ G n-p . Conditioning two times, first by H i,p-k and next by H i,p-+1 , and using the branching Markov property, we get

E Xi [M G p-( fn )M G p-k ( fn )] = 2 -k-1 E Xi [M G p-( fn )M G p-(g k, ,n )],
where we set g k, ,n = Q -k-1 P fn ⊕ Q -k-1 P fn . Next, conditioning by H i,p-and using the branching Markov property, we get

E Xi [M G p-( fn )M G p-(g k, ,n )] = E Xi [M G p-(P( fn g k, ,n ) -P fn Pg k, ,n )] + E Xi [M G p-(P fn )M G p-(Pg k, ,n )].
From the foregoing and using (52), ( 54) and (25), it follows that:

V 2 (n) = V 5 (n) + V 6 (n), where V 5 (n) = |G n-p | -1 M Gn-p (H 5,n ) and V 6 (n) = |G n-p | -1 M Gn-p (H 6,n ), with H 5,n = 0≤k< 2 -k Q p-(P( fn g k, ,n ))1 { ≤p} and H 6,n = 0≤k< r≥0 2 -k+r Q p--r-1 P(Q r P fn ⊗ sym Q r Pg k, ,n )1 {r+ <p} .
First, we treat the term V 6 (n). Note that we have

Q r Pg k, ,n = 2 Q r+ -k P fn .
We set

h (n)
k, ,r = 2 r-k+1 Q p-1-(r+ ) P(Q r P fn ⊗ sym Q r+ -k P fn ) and h k, ,r = 2 r-k+1 µ, P(Q r P fn ⊗ sym Q r+ -k P fn ) .

We consider the following sums:

H [n] 6 = 0≤k< r≥0 h k, ,r 1 {r+ <p} and A 6,n = H 6,n -H [n] 6 = 0≤k< r≥0 (h (n) k, ,r -h k, ,r )1 {r+ <p} .
Using Lemma 10.1, we have for all 0 ≤ k < and r ≥ 0:

|P(Q r P fn ⊗ sym Q r+ -k P fn )| ≤ Q r+ -k P fn ) ∞ P(Q r P fn ⊗ sym 1) ∞ ≤ Q r+ -k P fn ) ∞ Q r+1 P fn ∞ ≤ C h 3d n . (26)
Moreover, using (5) and Lemma 10.1, we have for all 0 ≤ k < and r ≥ 1:

|P(Q r P fn ⊗ sym Q r+ -k P fn )| ≤ C α 2r+ -k Q(Pf n ) 2 ∞ P(V ⊗ sym V ) ≤ C α 2r+ -k h 3d n P(V ⊗ sym V ). ( 27
)
Distinguishing the cases r = 0 and r ≥ 1 and using (26), ( 27), (iv) of Assumption 1.2 and (3), we get, for some g 1 , g ∈ F ,

|H 6,n -H [n] 6 | = 0≤k< |h (n) k, ,0 -h k, ,0 | 1 { <p} + 0≤k< r≥1 |h (n) k, ,r -h k, ,r | 1 {r+ <p} ≤ C 0≤k< 2 -k α p--1 P(P fn ⊗ sym Q -k P fn ) ∞ V + C h 3d n 0≤k< r≥1 2 r-k α 2r+ -k Q p--r-1 P(V ⊗ sym V ) 1 {r+ <p} ≤ C h 3d n (V + ( 0≤k< r≥1 2 r-k α 2r+ -k 1 {r+ <p} ) g 1 ) ≤ C h 3d n a n g, (28)
where

a n =      1 if 2α 2 < 1 p if 2α = 1 (2α 2 ) p if 2α 2 > 1.
Using (28), we find that 6), [START_REF] Parzen | On estimation of a probability density function and mode[END_REF] and that g ∈ L 1 (µ), we get

|V 6 (n) -H [n] 6 | ≤ |G n-p | -1 M Gn-p (|H 6 (n) -H [n] 6 |) ≤ C a n h 3d n |G n-p | -1 M Gn-p (g). Using (
lim n→+∞ |V 6 (n) -H [n]
6 | = 0 a.s. and in L 2 .

Next, as for (28), using (26), ( 27) and that F ⊂ L 1 (µ), we find that |H

(n) = (V 6 (n) -H [n] 6 ) + H [n] 6 , we conclude that lim n→∞ E[V 6 (n) 2 ] = 0.
Next, we treat the term V 5 (n). We have

V 5 (n) = (V 5 (n) -H [n] 5 ) + H [n]
5 , where

H [n] 5 = 0≤k< 2 -k µ, P( fn g k, ,n ) 1 { ≤p} .
Using Lemma 10.1, we get, for all 0 ≤ k < ,

|P( fn g k, ,n ) -µ, P( fn g k, ,n ) | ≤ C Q -k-1 P fn ∞ P fn ∞ ≤ C h d n 1 {k= -1} + C h 2d n 1 {k≤ -2} .
Using the latter inequality and distinguishing the cases = p and ≤ p -1, we find that

|V 5 (n) -H [n] 5 | ≤ C (2 -p h d n + h 2d n + p h 2d n ). This implies that lim n→+∞ |V 5 (n) -H [n]
5 | = 0 a.s. and in L 2 . Next, using Lemma 10.1, we have

|H [(n)] 5 | ≤ 0≤k< 2 -k | µ, P( fn g k, ,n ) |1 { ≤p} ≤ C >0 (2 -+1 h 2d n + -2 k=0 2 -k α -k h 3d/2 n ) ≤ C h 3d/2 n .
This implies that lim n→∞ H [(n)] 5

= 0 and then that lim n→∞ V 5 (n) = 0 a.s. and in L 2 .

Finally, since V 2 (n) = V 5 (n) + V 6 (n), it follows from the foregoing that lim n→∞ E[V 2 (n) 2 ] = 0 and this ends the proof.

Now we treat the term

V 1 (n). Recall V 1 (n) = |G n | -1 i∈Gn-p p k=0 E Xi [M G p-k ( fn ) 2 ].
We have the following convergence. Lemma 6.5. Under the assumptions of Theorem 3.3, we have the following convergence:

lim n→∞ V 1 (n) = 2 K 0 6 2 µ (x, x 0 , x 1 ) in probability.
Proof. Let k ∈ {0, . . . , p} and i ∈ G n-p . Conditioning by H i,p-k and using the branching Markov property, we get

E Xi [M G p-k ( fn ) 2 ] = E Xi [M G p-k (P( f 2 n ) -(P fn ) 2 )] + E Xi [(M G p-k (P fn )) 2 ]
. Using the latter inequality and the fact that P(

f 2 n ) -(P fn ) 2 = P(f 2 n ) -(Pf n ) 2 , we get V 1 (n) = V 3 (n) + V 4 (n) -V 7 (n), where V 3 (n) = |G n | -1 i∈Gn-p p k=0 E Xi [M G p-k (P(f 2 n ))]; V 7 (n) = |G n | -1 i∈Gn-p p k=0 E Xi [M G p-k ((Pf n ) 2 )]; V 4 (n) = |G n | -1 i∈Gn-p p k=0 E Xi [(M G p-k (P fn )) 2 ].
First we treat V 7 (n). We set

H 7,n = p k=0 2 -k Q p-k ((Pf n ) 2 ).
Using (52), we have

V 7 (n) = |G n-p | -1 M Gn-p (H
(n)| ≤ C (2 -p h d n + h 2d n ). It then follows that lim n→∞ V 7 (n) = 0 in probability. Next, we treat the term V 3 (n). We set A 3,n = H 3,n -H [n] 3 , with: (29) H 3,n = p k=0 2 -k Q p-k (P(f 2 n )) and H [n] 3 = p k=0 2 -k µ, P(f 2 n ) = 2(1 -2 -p-1 ) µ, P(f 2 n ) .
We set g n = P(f 2 n )µ, P(f 2 n ) . Using ( 5) and Lemma 10.1, we have

|V 3 (n) -H [n] 3 | ≤ |G n-p | -1 M Gn-p (2 -p |g n |) + M Gn-p ( p-1 k=0 2 -k |Q p-k-1 (Q(g n ))|) ≤ C 2 -p + C |G n-p | -1 M Gn-p ( p-1 k=0 2 -k α p-k Q(P(f 2 n )) ∞ V ) ≤ C 2 -p + C a n |G n-p | -1 M Gn-p (V ), (30)
where

a n =      2 -p if 2α < 1 p 2 -p if 2α = 1 α p if 2α > 1.
Using [START_REF] Comte | On a nadaraya-watson estimator with two bandwidths[END_REF] and the fact that lim n→+∞ a n = 0, we find that

lim n→+∞ C 2 -p + C a n |G n-p | -1 M Gn-p (V ) = 0 a.s. and in L 2 .
From (30), this implies that

lim n→+∞ |V 3 (n) -H [n] 3 | = 0 in probability.
Using Lemma 10.2, we get

lim n→∞ H [n] 3 = lim n→∞ 2 (1 -2 -p-1 ) µ, P(f 2 n ) = 2 K 0 6 2 µ (x, x 0 , x 1 ).
From the foregoing, we conclude that lim n→∞ V 3 (n) = 2 K 0 6 2 µ (x, x 0 , x 1 ) in probability.

Finally, we treat the term V 4 (n). Using (53), we have

V 4 (n) = V 8 (n) + V 9 (n), where V 8 (n) = |G n-p | -1 M Gn-p (H 8,n ) and V 9 (n) = |G n-p | -1 M Gn-p (H 9,n ), with H 8,n = p k=0 2 -k Q p-k ((P fn ) 2 ) and H 9,n = k≥0, ≥0 2 -k+ Q p-k--1 P(Q P fn ⊗ Q P fn )1 {k+ <p} .
Writing As a result, we find that

H [n] 8 = p k=0
lim n→∞ |V 8 (n)| = lim n→∞ |V 9 (n)| = 0 in probability. Since V 4 (n) = V 8 (n) + V 9 (n), we conclude that lim n→∞ V 4 (n) = 0 in probability. Finally, since V 1 (n) = V 3 (n) + V 4 (n) -V 7 (n)
, the result of the Lemma follows from the foregoing.

As a consequence of (24), Lemmas 6.3, 6.4 and 6.5, we have the following result.

Lemma 6.6. Under the assumptions of Theorem 3.3, we have the following convergence:

lim n→∞ V (n) = 2 K 0 6
2 µ (x, x 0 , x 1 ) in probability.

We now check the Lindeberg condition using a fourth moment condition. We set:

R 3 (n) = i∈Gn-p n E ∆ n,i (f n ) 4 .
Lemma 6.7. Under the assumptions of Theorem 3.3, we have that lim n→∞ R 3 (n) = 0

Proof. We have

R 3 (n) ≤ 16 (p + 1) 3 |G n | -2 i∈Gn-p p =0 E[(M iG p-( fn )) 4 ] ≤ 128 (p + 1) 3 |G n | -2 i∈Gn-p p =0 E[(M iG p-(f n -Pf n )) 4 ] (31) + 128 (p + 1) 3 |G n | -2 i∈Gn-p p =0 E[(M iG p-(P fn )) 4 ],
where we used that ( r k=0 a k ) 4 ≤ (r + 1) 3 r k=0 a 4 k for the two inequalities (resp. with r = 1 and r = p), Jensen inequality for the first inequality and the decomposition f n = (f n -Pf n ) + Pf n for the last inequality. For the second term of the right hand side of (31), we follow the proof of Lemma 5.6 in [START_REF] Bitseki Penda | Central limit theorem for kernel estimator of invariant density in bifurcating markov chains models[END_REF] and Lemma 4.7 in [START_REF] Bitseki Penda | Central limit theorem for bifurcating markov chains under pointwise ergodic conditions[END_REF] to find that

n 3 |G n | -2 i∈Gn-p p =0 E[(M iG p-(P fn )) 4 ] ≤ C n 5 (2 -n+p h n 1 {2α 2 ≤1} + 2 -n+p (2α 2 ) 2p h 6d n 1 {2α 2 >1} ),
and using [START_REF] Parzen | On estimation of a probability density function and mode[END_REF], this implies that (32) lim

n→+∞ n 3 |G n | -2 i∈Gn-p p =0 E[(M iG p-(P fn )) 4 ] = 0.
We are now going to treat the first term of (31). Since P(f n -P(f n )) = 0, we have, (see Remark 2.3 in [START_REF] Bitseki Penda | Deviation inequalities, moderate deviations and some limit theorems for bifurcating Markov chains with application[END_REF] for more details), Using (53), the fact that P((f n -Pf n ) 2 ) ≤ P(f 2 n ) for the first inequality and Lemma 10.1 for the second inequality, we get

ν, Q n-p h n, ≤ 2 p-ν, Q n-(P(f 2 n )) 2 + p--1 k=0 2 p-+k ν, Q n--k-1 P(Q k P(f 2 n ⊗ 2 ) ≤ C 2 p-+ C 2 2(p-) .
The latter inequality implies that (36)

n 3 2 -n-p p =0 ν, Q n-p h n, ≤ C n 3 (2 -n + 2 -n+p ).
From (34), ( 35) and (36), we conclude that lim n→∞ R 3,1 (n) = 0. Finally, from (31) and (32), this proves that lim n→∞ R 3 (n) = 0.

For A n = G n , we follows exactly the same lines and this ends the proof.

Lemma 8.2. Under Assumptions of Theorem 4.2, we have

(|A n |h 3d n ) 1/2 ( µ An (xx 0 x 1 ) µ An (x) - µ (xx 0 x 1 ) µ An (x) ) (d) ----→ n→∞ G,
where G is a centered Gaussian real-valued random variable with mean 0 and variance σ 2 = K 0 6 2 P(x, x 0 , x 1 )/µ(x).

Proof. This is a direct consequence of Theorem 3.5, (43) and Slutsky's Lemma.

9. Proof of (16)

We set f 0h (y) = h -1 K 0 (h -1 (xy)) and recall µ, f 0h = R f 0h (y)µ(y)dy. Using the decomposition

µ(x) -µ(x) = 1 |G n | u∈Gn f0h (X u ) + µ, f 0h -µ(x),
we obtain the following biais-variance type decomposition.

(45)

E ( µ(x) -µ(x)) 2 ≤ 2(|G n |) 2 E ( u∈Gn f0h (X u )) 2 + 2 µ, f 0h -µ(x) 2 .
Using (53) and the fact that the process L(X ∅ ) = µ (which implies that µQ = µ), we get (46)

|G n | -2 E ( u∈Gn f0h (X u )) 2 ≤ 2 µ, f 2 0h |G n | + 1 |G n | n-1 k=1 2 k µ, (Q k-1 (Q f0h )) 2 .
We now plan to use (15) with f = Qf 0h . For all y ∈ R, we get, after the change of variable t = h -1 (xz) and the use of the first-order Taylor's expansion, We then have that In this section, we recall useful results on BMC from Bitseki-Delmas [START_REF] Bitseki Penda | Central limit theorem for bifurcating markov chains[END_REF].

Lemma 10.3. Let f, g ∈ B(S), x ∈ S and n ≥ m ≥ 0. Assuming that all the quantities below are well defined, we have:

E x [M Gn (f )] = |G n | Q n f (x) = 2 n Q n f (x), (52)
E x M Gn (f ) 2 = 2 n Q n (f 2 )(x) + n-1 k=0 2 n+k Q n-k-1 P Q k f ⊗ Q k f (x), (53) E x [M Gn (f )M Gm (g)] = 2 n Q m gQ n-m f (x) (54) + m-1 k=0 2 n+k Q m-k-1 P Q k g ⊗ sym Q n-m+k f (x).

PPP

  estimation of P using the rule of thumb type method with n = 13; a0 = 0.7; a1 = 0.5;x0 x1The projection of PG n (left) and P (right) on the plane x = 3. x0 x1 P estimation of P using the rule of thumb type method with n = 15; a0 = 0.7; a1 = 0.5; x0 x1The projection of PG n (left) and P (right) on the plane x = 3.

Figure 1 x0 x1 PPP

 1x1 Figure 1

Figure 2 6.

 2 Figure 2

PPP

  estimation of P using the rule of thumb type method with n = 13; a0 = 1.2; a1 = 0.7;x0 x1The projection of PG n (left) and P (right) on the plane x = 3. x0 x1 P estimation of P using the rule of thumb type method with n = 15; a0 = 1.2; a1 = 0.7;x0 x1The projection of PG n (left) and P (right) on the plane x = 3.

Figure 3 x0 x1 PPP

 3x1 Figure 3

  Figure 4

  Figure 6

  projection of PGn (black) and P (green) on the line x = 3; x1 = 0. We use the cross validation method and n = 13; a0 = 1.2; a1 = 0.7; projection of PGn (black) and P (green) on the line x = 3; x1 = 0.We use the cross validation method and n = 15; a0 = 1.2; a1 = 0.7; x0 P (black) and P (green)

  C a n h 3d n . Using[START_REF] Parzen | On estimation of a probability density function and mode[END_REF], we get lim n→+∞ H[n] 6 = 0. Now, since we can write V 6

  (33) E x [(M G p-(f n -Pf n ))4 ] ≤ g n, (x) + 6 h n, (x), with:g n, (x) = E x [M G p-(P((f n -Pf n ) 4 ))] and h n, (x) = E x [(M G p-(P((f n -Pf n ) 2 ))) 2 ]. We set R 3,1 (n) = (p + 1) 3 2 -2n i∈Gn-p p =0 E[(M iG p-(f n -Pf n )) 4 ].Using the branching Markov property and (33) for the first inequality and (52) for equality, we getR 3,1 (n) ≤ C n 3 2 -2n p =0 E[M Gn-p (g n, )] + C n 3 2 -2n p =0 E[M Gn-p (h n, )] = C n 3 2 -n-p p =0 ν, Q n-p g n, + C n 3 2 -n-p p =0 ν, Q n-p h n, .(34)Using Lemma 10.1, we getν, Q n-p g n, = 2 p-ν, Q n-(P((f n -Pf n ) 4 )) ≤ 2 p-QP((f n -Pf n ) 4 ) ∞ ≤ C h -3d n 2 p-.The latter inequality implies that (35)n 3 2 -n-p p =0 ν, Q n-p g n, ≤ C n 3 (2 n h 3d n ) -1 .

(

  Qf 0h ) (y) = 1 h R K 0 (h -1 (xz)) ∂Q ∂y (y, z)dz = R K 0 (z) ∂Q ∂y (y, xhz) dz = ∂Q ∂y (y, x) R K 0 (z) dz + O(1) = ∂Q ∂y (y, x) + O(1).

(√ 2πσ 2 e - 1 / 2 2 1 ( 1

 21221 Qf 0h ) ∞ = sup y∈R ∂Q ∂y (•, x) + O(1) = 1 +O(1).Using the latter equality and (15), we get, for all k ≥ 1,(47) |Q k-1 (Q f0h )|(y) ≤ 1 √ 2πσ 2 e 1/2 (1a) (σ(1 + a) -1 + |y|)a k-1 + O(1).Recall µ is the Gaussian law N(0, σ 2 (1a 2 )). Using (47) and (σ(1+ a) -1 + |y|) 2 ≤ 2 σ 2 (1 + a) -2 + 2 y 2 , we get (48) µ, (Q k-1 (Q f0h )) 2 ≤ 1 eπσa) 2 + 1 1a 2 a 2(k-1) + O(1) ≤ C a,σ a 2 a 2k + O(1).

Lemma 10 . 2 .

 102 Let (h n , n ∈ N) be a sequence of positive numbers converging to 0 as n goes to infinity. Let g : R d → R be a measurable function such thatR d |g(x)|dx < +∞. Let f : R d → R be a measurable function such that f ∞ < +∞, R d |f (y)| dy < +∞ and lim |x|→+∞ |x|f (x) = 0. Define g n (x) = h -d n R d f (h -1 n (x -y))g(y)dy. Then, we have at every point x of continuity of g, lim n→+∞ g n (x) = g(x) R f (y)dy.

  7,n ). Using Lemma 10.1 and distinguishing the cases k = p and k ≤ p -1, we get |V 7

We can now use Theorem 3.2 and Corollary 3.1, p. 58, and the remark p. 59 from [START_REF] Hall | Martingale limit theory and its application[END_REF] to deduce from Lemmas 6.6 and 6.7 that ∆ n (f n ) converges in distribution towards a Gaussian real-valued random variable with deterministic variance σ 2 . Using (20) and Lemmas 6.1 and 6.2, we then deduce Theorem 3.3 for A n = T n .

For A n = G n , we have

and for all i ∈ G n-p ,

Following exactly the proof of Lemma 6.2, 6.6 and 6.7, we get the result for this case. We note that for A n = G n , the factor 2 is missing in the asymptotic variance. This come from the fact that here, H

[n]

3 defined in (29) is simply equal to µ, P(f 2 n ) .

7. Proof of Theorem 3.5

We begin the proof with A n = T n . From [START_REF] Hall | Martingale limit theory and its application[END_REF], we have

where the bias term B n (xx 0 x 1 ) is defined by

Since lim n→∞ (|G n |/|T n |) 1/2 = 1/ √ 2, from Theorem 3.3, it suffices, to obtain the result of Theorem 3.5, to prove that lim n→∞ B n (xx 0 x 1 ) = 0. Using the Taylor expansion and Assumption 3.4, one can prove that (see [START_REF] Bitseki Penda | Central limit theorem for bifurcating markov chains[END_REF] for more details)

Since lim n→∞ |T n |h 2s+3d n = 0, we conclude that lim n→∞ B n (xx 0 x 1 ) = 0 and this ends the proof for A n = T n .

For A n = G n the proof follows exactly the same lines.

Proof of Theorem 4.2

First of all, we have the following decomposition:

Then, the proof of Theorem 4.2 is a direct consequence of the previous decomposition and Lemmas 8.1 and 8.2 below.

Lemma 8.1. Under Assumptions of Theorem 4.2, we have

Proof. We consider the function g n defined on S by: g n (y) = h

for all y ∈ S. We begin the proof with A n = T n . We set gn = g nµ, g n . We have the following decomposition:

Using the fact that K 0 is bounded, integration by parts and Assumption 1.8, we have the following upper bounds:

Next, from Minkowski's inequality, we have

Using (53), ( 5) and (38), we get:

The latter inequality implies that

Using (39), ( 40) and ( 11), we deduce that

Next, using Taylor expansion and Assumption 3.4, we get (see [START_REF] Bitseki Penda | Central limit theorem for bifurcating markov chains[END_REF] for more details) 37), ( 41) and (42), we deduce that (43) lim n→∞ µ Tn (x) = µ(x) in probability.

We further deduce that

Using [START_REF] Parzen | On estimation of a probability density function and mode[END_REF], the latter inequality implies that 43), (44) and using Slutsky's Lemma, we get

Now, (48) and (46) implies that

Putting the latter inequality into (45), we obtain

Finally, it is very standard to get asymptotic equivalence of the second and the third term of the right hand side of (49) (see for e.g. [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF], Section 3.3.1 for more details). This ends the proof of (16).

Appendix

First, we give some useful upper bounds. We recall that S = R d . Recall f n defined in [START_REF] Gyori | Hypothesis testing for markov chain monte carlo[END_REF].

Lemma 10.1. Under Assumption (1.8), we have:

Proof. Using a change of variables, we have, for all y ∈ S:

|K 0 |(y 0 ) |K 0 |(y 1 )P(y, x 0h y 0 , x 1h y 1 )dy 0 dy 1 .

This implies that

n . From (50) and using again a change of variables, we have, for all t ∈ S:

|K 0 |(y) |K 0 |(y 0 ) |K 0 |(y 1 ) P(xh y, x 0h y 0 , x 1h y 1 ) Q(t, xh y)dy dy 0 dy 1 .

This implies that (51)

As for (51), we have

n . Now, following the same ideas, we easily get the others upper bounds.

We recall the following result due to Bochner (see [START_REF] Parzen | On estimation of a probability density function and mode[END_REF]Theorem 1A] which can be easily extended to any dimension d ≥ 1).