Loann D Desboulets
email: loann.desboulets@univ-amu.fr

Non-parametric Variable Selection on Non-linear Manifolds

Keywords: Non-parametric, Non-linear Manifolds, Variable Selection, Neural Networks 1

In this paper, I investigate a new non-parametric variable selection framework.

To extend the usual non-parametric model, I consider non-linear manifolds which are more flexible. Non-linear manifolds are represented by function compositions, allowing more complex dependences in the data. Based on two manifold approximation techniques, k-nearest neighbours and auto-encoder neural networks, I propose two different procedures to perform non-parametric variable selection. The two methods are complementary, the former being a local estimator, while the latter is a global estimator.

Introduction

Non-parametric models are popular modelling techniques, that allows a great flexibility.

When the number of explanatory variables is large, variable selection come in handy. Variable selection is meant to infer a subset of the most pertinent explanatory variables from the nonparametric model y i = m(x i ) + i ∀i ∈ 1, .., n.

(

) 1
Where y is the dependent variable, x is a vector of covariates, is a noise term and m is an unknown function. Selection is useful in reducing the dimension of the estimator m(•), which is a serious issue for high dimensional non-parametric models, known as the curse of dimensionality. There exists plenty of methods to perform non-parametric variable selection, see [START_REF] Mehmood | A review of variable selection methods in partial least squares regression[END_REF]; [START_REF] Desboulets | A review on variable selection in regression analysis[END_REF] for recent reviews. However, there might exists more complex correlations, that cannot be modelled directly using functions.

And in equation 1, the unknown model m(•) has to be a function. The Swiss Roll data of [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] is an instance of such complex correlations. To get a more general approach, I consider non-linear manifolds which are more flexible than functions.

A manifold is often represented using a function composition. This can represent more complex relationships than a single function. Also in the manifold equation, there is also no distinction between explanatory variables and outcome variables. Both are embed into a single matrix (X 1 , X 2 , ..., X p ) ≡ X ∈ R n×p . This can be useful, in situations where there is a priori no obvious reason to choose a particular variables as being more dependent than the others. It can occurs in many areas of Economics, especially in Finance, for example when regressing some returns on other returns. The manifold M is written as follow

M :={x ∈ R p |x = f • g(x)} g : R p → R r , f : R r → R p , p > r (2)
The projection function g maps the p-dimensional data onto the lower-dimensional coordinates of the r-manifold. Then, the immersion function f is mapping back to the original space. The data are assumed to be i.i.d. and stationary. The manifold selection problem is to find the variables which have the greatest correlations with all the other variables, expressed in terms of marginal contributions. Obviously, I have to discard the correlation between a variable and itself. Its projection through the manifold cannot be the identity function. The "true" set of these correlated variables is denoted S * and defined as follow

S * := {j : ∃j = j, ∂X j ∂X j = 0} ⇔ {j : ∂(f • g) ∂X j = 1} |S * | ≤ p (3)
The contribution of this paper is to propose two manifold selection operators to infer the set of "active variables" S * , from equation 3. The first is a locally linear estimator, based on the k-nearest neighbours algorithm named LLMS, for Local Linear Manifold Selection.

The second is a global estimator, based on auto-encoders neural networks [START_REF] Kramer | Nonlinear principal component analysis using autoassociative neural networks[END_REF], named DAMS for Diagonal Auto-encoder Manifold Selection. The two approaches rely on completely different mechanisms but solve almost the same problem. The local estimator does not rely on continuity everywhere, contrary to the global one. This means f , and or g, can be only locally continuous. Moreover, the local approach has low computational complexity but might be sensitive to noisy data. The global one is exactly the opposite. These reasons motivate for having two distinct approaches.

In general, any selection operator relies on a regression estimator, that approximates an associated model. Thus, selection on manifold has to rely on an estimator that approximates the manifold. A large literature is devoted to inference of manifold in the data, referred to as principal manifolds or unsupervised manifold learning. They use non-parametric estimators such as for instance kernels [START_REF] Hastie | Principal curves[END_REF]; [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]; [START_REF] Mika | Kernel pca and de-noising in feature spaces[END_REF], or neural networks [START_REF] Kohonen | The self-organizing map[END_REF]; [START_REF] Vincent | Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion[END_REF], and nearest neighbours [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF]; [START_REF] Saul | Think globally, fit locally: unsupervised learning of low dimensional manifolds[END_REF]. A useful introduction to manifold learning can be found in [START_REF] Gorban | Principal graphs and manifolds[END_REF].

Obviously these are regression techniques, they are not doing variable selection per se. They only estimate the manifold. The literature on manifold selection is quite recent, but vary across topics. On the one hand, there exists several papers for classification problems and image recognition [START_REF] Ding-Cheng | Detecting local manifold structure for unsupervised feature selection[END_REF][START_REF] Li | Manifold regularized multi-view feature selection for social image annotation[END_REF][START_REF] Tang | Unsupervised feature selection via latent representation learning and manifold regularization[END_REF][START_REF] Doquet | Agnostic feature selection[END_REF]. Typically, they assume the explanatory variables lie on a manifold. This manifold represents the only source of variability in the data. It summarizes all the information contained in the data, and so it reduces the dimension of the estimator, avoiding de facto the curse of dimensionality. It is a non-parametric equivalent to Principal Component Regression [START_REF] Jolliffe | A note on the use of principal components in regression[END_REF]. Once extracted, using the manifold instead of the raw data produces a more parsimonious classifier, usually with better performance. The paper of [START_REF] Ding-Cheng | Detecting local manifold structure for unsupervised feature selection[END_REF] is actually very close to the first of the two methods I propose (LLMS). The main difference is they use a score based on reconstruction error to assess the importance of a variable. They focus on finding which variables is correlated to the binary outcome, which is not what I do here. I try to find which variables are correlated with all the others, not a single one. On the other hand, there is almost no literature on manifold selection for regression purposes, i.e. only with continuous data. The only method that is close to the one I propose, is the Nonparanormal (hereafter NPN) of [START_REF] Lafferty | Sparse nonparametric graphical models[END_REF]. In their paper, they assume non-parametric transformations of the variables are distributed as a multivariate Gaussian. Once the transformations are estimated, the selection is performed using the Graphical LASSO Friedman et al. (2008). However, their method is restricted to monotonically increasing transforms. This is a particular kind of manifolds. Hence, their method is not as flexible as the method I propose here.

The manifold model is a general and flexible statistical model. That is desirable, because assuming an inappropriate statistical model might not lead to consistent selection. Intuitively, if a selection method assumes that the data is normally distributed, while it is not, it is unlikely the estimated set of correlated variables is the true one. More formally, let us define a selection operator as S = σ(X, P ). This operator depends on a statistical model, described by a sample X, and a set of parametrized probability distributions P = {P ω : ω ∈ Ω}.

For instance, the famous LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] is based on a linear statistical model

{Y = Xβ : β ∈ R p }.
For the sake of simplicity, hereafter the notation for the distribution is shortened to its parameters. The notation for the LASSO is therefore σ LASSO (X, β). This LASSO achieves sparsity by solving min β Y -Xβ 2 2 + θ β 1 , where θ is a hyper-parameter that determines the amount of regularisation. Generally speaking, the Manifold Selection (MS) estimator is denoted as σ M S (X, f • g). So, using simulations, I show there must be some correspondence between the "true" data generating process and the assumed statistical model, otherwise the selected set is not consistent. Only the two methods I propose infer S * correctly when the data is distributed on a manifold. Without further restrictions, the identification of S * is not possible. Also, the function composition is not unique. The condition on the domains of the functions p > r reduces the function space for sure, but is not sufficient to obtain uniqueness. To identify S * , what is required is the identification of the partial derivatives, not f

• g itself. Still, if f • g is unique,
so are its derivatives. The fundamental problem here, is that it is always possible to find a function, as complex as needed, that links two variables. More generally, a smooth curve (r = 1) can pass through every points in the data, provided it is non-linear enough. One simple and usual solution for non-parametric models, is to impose a regularity condition on the total variation of f • g. In the literature, a penalty on the second order derivatives is often used. Indeed, this limits, by definition, the variability of the function. But there exists other ways to prevent from the aforementioned problem. In the case of the global estimator, based on neural networks, it suffices to have "not too large" intermediate layers.

If the activation functions are monotonic, each supplementary neurons in the intermediate layer brings another possible extrema to the fitted function. Limiting the number of neurons therefore limits the variability. For the local estimator, the solution is the local linear models themselves, as a linear function has no variability. Notice that this is true as long as the parameter k is large enough: k ≥ p + 1. Otherwise there no identification of the linear models.

Another identification problem arises if we consider the variables to be measured with noise, which is often the case empirically.

M :={x ∈ R p |x = f • g(x * ) + } g : R p → R r , f : R r → R p , p > r (4)
Here each variable is contaminated with a specific error term, therefore the dimension of the noise is the same as the dimension of the data. Another possibility is for the coordinates of the manifold to be noisy, but I do not take this one in consideration. In the general case, the conditions that ensures identification of the manifold are not known. A necessary but not sufficient condition is that true values of the data are independent from the noise x * ⊥ , and the noises are uncorrelated i.e. V[ ] is diagonal. There is no need to impose the same variance for each noise, nor this variance to be constant over all the distribution of the data.

The rest of the paper is organised as follow. In section 2, I develop the LLMS estimator, based on the k-nearest neighbours algorithm. Its theoretical properties are studied, and three different versions are proposed to improve selection. In section 3, I develop the DAMS estimator, based on neural networks. I study a new network structure to perform selection and extend it to auto-encoder networks. In section 4, I run simulations to compare existing methods, LLMS and DAMS. I show that both proposed methods can select variables consistently.

Local Linear Approach via k-Nearest Neighbours

Motivations for a Locally Linear estimator

The very definition of a manifold is that it is locally a Euclidean space, i.e. in a neighbourhood N of a given point x 0 defined in a small enough euclidean ball B x 0 of radius ν.

B x 0 ={X ∈ R p , X -x 0 2 2 ≤ ν} ⊂ R p N x 0 ={x ∈ B x 0 |x = f • g(x)}.
(5)

It is also worth noticing that the union of all the possible neighbourhoods is the entire manifold

M = ∀x 0 N x 0 . (6)
The data matrix X is a sample from M, and the neighbourhoods are sub-samples of X. In a neighbourhood of x 0 , a linear approximation of the manifold is sufficient as the error term is limited by definition to be a small value:

O((X -x 0 ) 2 ) ≈ X -x 0 2 2 ≤ ν.
Following a first order Taylor expansion of the manifold:

f • g(x) = f • g(x 0 ) + ∇ f •g (x 0 )(X -x 0 ) + O((X -x 0 ) 2 ), (7)
the definition of a neighbourhood can be changed to

N x 0 = {x ∈ B x 0 |x ≈ f • g(x 0 ) + ∇ f •g (x 0 )(X -x 0 )} ⇔N x 0 = {x ∈ B x 0 |x ≈ µ x 0 + xβ x 0 } (8)
Following equations 6 and 8, X is the union of all the neighbourhoods and the function composition f •g is approximated by the union of all the local linear models. So the manifold selection can be rewritten as follow:

σ(X, f • g) = σ( ∀x 0 N x 0 , ∀x 0 µ x 0 + xβ x 0 ). (9)
Now, if I assume separability of the selection operators, I can write the manifold selection as a union of local linear selection:

σ( ∀x 0 N x 0 , ∀x 0 µ x 0 + xβ x 0 ) = ∀x 0 σ(N x 0 , µ x 0 + xβ x 0 ). ( 10
)
To illustrate this local separability assumption, let's take a practical example. Assume we have three variables, x, y, z, such that x and y lie on a circle and z is a noise variable. We can write the associated manifold using a function composition. The projection function g maps the two dimensional coordinates x and y onto the one-dimensional circle, and keeps z as is. The immersion function f is basically the inverse of g.

g : R 3 → R 2 (x, y, z) → (t := atan2(x, y), z). f : R 2 → R 3 (t, z) → (sin(t), cos(t), z). ( 11
)
with atan2 : R 2 → R 1 (x, y) → atan(y/x)1 {x>0} + π1 {x < 0 & y ≥ 0} -π1 {x < 0 & y < 0} + π 2 1 {x = 0 & y > 0} - π 2 1 {x = 0 & y < 0} (12)
In matrix form, we have:

X :=(x, y, z), X =f • g(X). (13)
The manifold, which is a cylinder, formed by the three variables x, y, z is depicted in Figure 1. The "true set" of variable S * := {j : ∃j = j, ∂X j ∂X j = 0} requires to compute partial derivatives of each variable w.r.t to all the others. These derivatives are easier to calculate using the implicit form of the manifold equation. Let's define the low-coordinates of the manifold as follow:

t ∼ U[-π, π], ∼ U[-1, 1]. ( 14
)
Then the manifold can be written in implicit form as follow:

x = sin(t), y = cos(t), z = .

(15) Fig. 1. Example of a non-linear manifold formed with three variables x, y, z

These equations can also be written in terms of the observed variables x, y, z only:

x = sin(cos -1 (y)) = (1 -y 2 ), y = cos(sin -1 (x)) = (1 -x 2 ), z = (16)
These formulations are equivalent to the formulation in equation 13. Now, we can work out the partial derivatives:

∂x ∂y = - y (1 -y 2 ) =          = 0 if y = {0, -1, +1} = 0 if y = 0 undefined if y = {-1, +1} ∂y ∂x = - x (1 -x 2 ) =          = 0 if x = {0, -1, +1} = 0 if x = 0 undefined if x = {-1, +1} ∂x ∂z = ∂z ∂x = ∂y ∂z = ∂z ∂y = 0 (17)
Or in matrix form:

∂X ∂X =    1 -y √ (1-y 2 ) 0 -x √ (1-x 2 ) 1 0 0 0 1    (18)
The true set of variable is therefore:

⇒ S * := {j : ∃j = j, ∂X j ∂X j = 0} = {1, 2} (19)
It is worth noticing here, what a selection operator assuming a linear model in the data would do. The plane that best fit the data in Figure 1, is flat in all variables. The correlation matrix of these data is a diagonal matrix, and so are its eigenvectors. So whatever the mechanism that drives the selection, e.g. penalty on the parameters, it would not select anything. This happens because the supports of all three variables are symmetric ([-1, +1]).

However, with a suitable estimate of the partial derivatives, we can define a selection estimator such that σ(X, f • g) = S * . However, an estimate of the partial derivatives on the whole support is hard to obtain. We can recover the true set, only using local estimates of the derivatives. For the sake of simplicity, let us choose 4 points t 0 and compute the local linear selection at these points. The whole process is illustrated in Figure 2.

t 0 ∈ π 4 , -π 4 , 3π 4 , -3π 4
x 0 = (sin(t 0 ), cos(t 0 ), z(t 0 )) (20)
The local derivatives are obtained by taking the first order Taylor expansion f (t) = f (t 0 ) + f (t 0 )(t -t 0 ) of the circle and then differentiating w.r.t to the other variables.

• t 0 = π 4 x = sin(t) ≈ + 1 √ 2 + 1 √ 2 t - π 4 y = cos(t) ≈ + 1 √ 2 - 1 √ 2 t - π 4 (21)
One can see the two approximations have the same coefficients, and only differ by the leading coefficient sign, therefore the derivatives between x and y are trivial. The derivatives involving z do not change, the variable z is not present in the approximations formulas. ∂x ∂y = -1,

∂y ∂x = -1, ∂x ∂z = ∂z ∂x = ∂y ∂z = ∂z ∂y = 0. ( 22
)
Thus, the selection set at this point is σ(N x 0 ( π 4 ) , µ x 0 ( π 4 ) + xβ x 0 ( π 4 ) ) = {1, 2}. The same logic is applied to the other points.

• t 0 = -π 4 x = sin(t) ≈ + 1 √ 2 + 1 √ 2 t + π 4 y = cos(t) ≈ + 1 √ 2 + 1 √ 2 t + π 4 → ∂x ∂y = -1 → ∂y ∂x = +1 → ∂x ∂z = ∂z ∂x = ∂y ∂z = ∂z ∂y = 0 σ(N x 0 ( -π 4 ) , µ x 0 ( -π 4 ) + xβ x 0 ( -π 4 ) ) = {1, 2} (23)
• t 0 = 3π x = sin(t) ≈ + 1 √ 2 - 1 √ 2 t - π 4 y = cos(t) ≈ + 1 √ 2 - 1 √ 2 t - π 4 → ∂x ∂y = +1 → ∂y ∂x = -1 → ∂x ∂z = ∂z ∂x = ∂y ∂z = ∂z ∂y = 0 σ(N x 0 ( 3π 4 ) , µ x 0 ( 3π 4 ) + xβ x 0 ( 3π 4 ) ) = {1, 2} (24)
• t 0 = -3π

4 x = sin(t) ≈ + 1 √ 2 - 1 √ 2 t + π 4 y = cos(t) ≈ + 1 √ 2 + 1 √ 2 t + π 4 → ∂x ∂y = +1 → ∂y ∂x = +1 → ∂x ∂z = ∂z ∂x = ∂y ∂z = ∂z ∂y = 0 σ(N x 0 ( -3π 4 ) , µ x 0 ( -3π 4 ) + xβ x 0 ( -3π 4 ) ) = {1, 2} (25)
In the end, the union of the local selection leads to the same set of selected variables: This is required for identification of the parameters. The low dimensional projection can be mapped back using the transposed of the same vectors. Doing so, produces the best linear approximation in a total least square sense.

σ(X, f • g) = ∀x 0 σ(N x 0 , µ x 0 + xβ x 0 ) = {{1, 2}, {1, 2}, {1, 2}, {1, 2}} = {1, 2} = S * . ( 26
x

≈ (x -µ x )V 1:r V 1:r + µ x . (27)
The last step is to define a selection operator based on these parameters. There have been several attempts in the literature to perform variable selection in that case. The Sparse PCA of [START_REF] Zou | Sparse principal component analysis[END_REF] uses a composite L1-norm and L2-norm constraint on the eigenvectors.

However it would requires the knowledge of r, which is not the case here. The Graphical LASSO of [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] is the natural extension of the LASSO selection operator to covariance matrices. It is a serious candidate to perform the local linear selection, but it does not operate on eigenvectors, it operates directly on the covariance matrix. Thus, it does not take the dimension reduction into accounts. Also, even if there have been a very great development of algorithms to solve the graphical LASSO efficiently, the computational burden of many such estimations is large. Finally, the Penalized Matrix Decomposition (PMD) [START_REF] Witten | A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis[END_REF], which uses only L1-norm penalty on the eigenvectors, but also requires the knowledge of r. All the mentioned methods uses norm penalties to achieve sparsity, which have no analytical solutions and therefore are longer to compute.

Generally speaking, selection directly on all the eigenvectors is not viable because they are normalized. This means that elements of V r+1:p are on the same scale as V 1:r , between 0 and 1. And that affects the selection. Remind that r is not known nor estimated. Nevertheless, the eigenvalues L (in diagonal matrix form) do carry some information about r. They measure the lengths of the principal directions of a plane. Locally the r-manifold is a rdimensional plane embedded into R p . Eigenvalues are sorted in decreasing order, so there should be a significant difference between L r,r and L r+1,r+1 . So a way to account for the uncertainty on r is to use the eigenvectors weighted by the square root of eigenvalues: VL 1/2 .

There have also been propositions in the literature to try to infer r. In which case, I could use that estimate and keep only the corresponding columns of V. Still, the solution with VL 1/2 presented before is simpler and faster to compute. Then, to achieve sparsity I will only use a simple threshold on each element of VL 1/2 . This is solution is possible, instead of using L1 or L2 norms, because both orthogonality and shrinkage are already there from the eigenvalue decomposition. Having normalized the eigenvectors is equivalent to a shrinkage, in the same way as a ridge estimator. Also, the eigenvectors are orthogonal by construction.

In a classical regression problem, the non-orthogonality of the parameters is the reason why the LASSO has no analytical solutions. Thus, I propose a new linear selection operator: the EigenThresholding (ET).

Σ = VLV V ij = V ij 1 |V ij L 1/2 j,j |>θ . (28)
The parameter θ is a penalty that controls the amount of sparsity. To the best of my knowledge, this selection estimator has not been proposed yet. Its properties as a linear selection estimator are examined in the simulations, along with the other linear selection operators. The standard approach, if I were to follow the literature, would be to use a regression operator to obtain an estimate of the manifold, and then derive a selection Ŝ. But the crucial point is that I perform selection without ever explicitly estimate the manifold. The full procedure is described in Algorithm 1.

Local Linear Manifold Selection

A neighbourhood around a given observation x i (a point on the manifold) contains k neighbours, and is defined as:

N i,k (x i ) = x ∈ X : x -x i ≤ x -x i (k) (29)
The hyper-parameter k needs to be tuned. In each of these neighbourhoods, I apply a selection operator σ(X, Ω). In the previous section, I have motivated a choice for a specific estimator. As it was discussed there is not a single estimator suitable for this task. So here I present the general method, using any linear selection operator. Most of them solve a constrained optimisation problem, using a penalty function on the parameters to be estimated.

The constraint ensures sparsity in the parameters, which induces the selection. The penalty function usually involves a penalty parameter θ, it controls the amount of sparsity.

I choose to decompose the correlation matrix for three reasons. First, it standardizes the coefficients associated with each variable. There is no bias towards a candidate with a larger scale. This would have been the case using a covariance matrix. Second, it makes no arbitrary assumption on the dependent variable. All candidates possibly affect each other. This would not be true using the simple linear regression. The optimization problem in each neighbourhood N i,k is:

R N i,k = Cor (N i,k ) argmin RN i,k R N i,k -RN i,k + penalty RN i,k , θ (30)
Then I record the selected variables:

ŜN i,k ,θ = j| RN i,k (j, s) = 0, ∀s = j .
The selections ŜN i,k ,θ are turned into boolean vectors:

Ψ N i,k ,θ = 1 j∈ ŜN i,k ,θ
There is a vector Ψ for each neighbourhood and for each value of θ.

Finally, I average the n boolean selections Ψ N i,k ,θ .

P θ = 1 n n i=1 Ψ N i,k ,θ .
The vector P θ measures for each candidate variable, the proportion of neighbourhoods in which it was selected. It is simple to see that P θ ∈ [0, 1] p . So I might see it as a "pseudoprobability", for a given variable, to belonging to the manifold.

We can connect the matrix P to very usual objects known in the variable selection literature as "selection paths". Usually, one represents the behaviour of the parameter of interest as a function of the penalty parameters θ. Here the same is made but with P, in the same spirit as the "probability paths" of StabSel [START_REF] Meinshausen | Stability selection[END_REF]). The exact definition of θ is narrowed to the selection estimator σ(X, Ω). For some estimator, there have been procedure proposed to choose θ optimally. This issue is discussed in section 2.4.

The algorithm for the LLMS procedure is presented in Algorithm 1. Algorithm 1 Local Linear Manifold Selection (LLMS1) input X, k, Θ, σ(X, Ω)

for i = 1 to n do Search for the neighbourhood N i,k (X i ) for θ ∈ Θ do Estimate the sparse correlation matrix RN i,k using σ(N i,k , Ω) Record active variables in ŜN i,k ,θ = j| RN i,k (j, s) = 0, ∀s = j . Turn Ŝ into boolean vectors Ψ N i,k ,θ = 1 j∈ ŜN i,k ,θ
end for end for

Compute the final selection estimator

P θ = 1 n n i=1 Ψ N i,k ,θ . output P θ

Randomization of Distances for Large Candidate Set

There is an issue when constructing neighbourhoods as the candidate set increases. I cannot measure distances onto the manifold directly because S * is unknown. I only observe distances in the R p space of the candidate variables. So, I can only compute noisy versions of the true distances. It disturbs the ordering of the distances and therefore the neighbourhood is not the same. The local linear assumption might not hold anymore, and as a consequence, it may affect the selection.

A solution to obtain measures of distances closer to the desired ones is to use a modified version of kNN known as Random-kNN (RkNN) [START_REF] Ho | Nearest neighbors in random subspaces[END_REF]. Here, I use random subsets

B ⊂ A ≡ {j ∈ N, j ≤ p} of the candidates variable set A for computing distances. D kN N = x -x i ∀x ∈ X A . N kN N i,k = x ∈ X : D kN N ≤ D kN N (k) D RkN N = x -x i ∀x ∈ X B . N RkN N i,k = x ∈ X : D RkN N ≤ D RkN N (k)
where the subscript (k) indicates the k th order statistics. max

N kN N S * -x i ≥ max N RkN N S * -x i
Proof. Let X be a matrix of p random variables with n observations, so that

x i = f • g(x i
) ∀i ∈ {1, ..., n}, x i ∈ X. Now let's divide the full set of candidates X into 3 nonoverlapping matrices X , Y and Z composed of respectively p x , p y and p z random vectors such p x + p y + p z = p. Only X ∈ S * , while (Y, Z) / ∈ S * ). Let x 0 ∈ X be a given random point on the manifold. Wlog, let x 1 ∈ X and x 2 ∈ X be respectively the first and second closest point to x 0 :

px (x 1 -x 0 ) 2 < px (x 2 -x 0 ) 2 < px (x i -x 0 ) 2 ∀i = {0, 1, 2}, (31)
where the power of vector is defined element-wise. W.l.o.g. we can write:

px (x 2 -x 0 ) 2 = px (x 1 -x 0 ) 2 + x with x > 0.
Random-kNN consists in computing distances from random subset taken from X, e.g. only X . Now let's define y ∈ Y and z ∈ Z, given points which are not on the manifold. Now lets take two sets of indices A and B such that indices of X that correspond to (X , Y) are in B and A contains all X, so B ⊂ A. Now let's compare the distances w.r.t to (x 0 , y 0 , z 0 ).

The condition so the ordering of distances in equation 31 is preserved is

px (x 1 -x 0 ) 2 + py (y 1 -y 0 ) 2 < px (x 2 -x 0 ) 2 + py (y 2 -y 0 ) 2 + pz (z 2 -z 0 ) 2
Which can be rewritten as

px (x 1 -x 0 ) 2 - px (x 2 -x 0 ) 2 < py (y 2 -y 0 ) 2 - py (y 1 -y 0 ) 2 + pz (z 2 -z 0 ) 2 . ( 32
)
The first quantity tends to zero as the manifold is more and more sampled (n → ∞).

lim n→∞ px (x 1 -x 0 ) 2 - px (x 2 -x 0 ) 2 ≡ x → 0.
The distances between two points from Y depends on its distribution. If we assume each vector of Y to be i.i.d. and having finite first and second moments, then their distances are also i.i.d. with finite first and second moment. We can then write

E py (y 2 -y 0 ) 2 = E py (y 1 -y 0 ) 2 = µ y
Then, it follows that lim py→∞ py

(y 2 -y 0 ) 2 - py (y 1 -y 0 ) 2 ≡ y → 0.
Finally, since the sum of distances among Z only grows with the number of elements it contains lim pz→∞ pz

(z 2 -z 0 ) 2 ≡ z → ∞.
So we can rewrite equation 32

xy < z .

The probability that this condition is true tends to one as the number of irrelevant variables grows lim py→∞ pz→∞

P ( z > x -y ) = 1.
But how can I ensure that S * ∈ B and m ≡ |B| < |A| = p ? This is an ill-posed problem as it would require the knowledge of S * .

One solution is to start by sampling the subset B from {j ∈ N, j ≤ p}, using a uniform distribution π. It makes no a priori on which variables are more likely to belong to the manifold. Then, I can update that sampling probability π, using selection probabilities from previous steps. As the process goes on, the subset will be drawn from most probably active variables. The ordering of the data may have an impact here, if the selection is biased in the first neighbourhoods. Still, the sampling probability might correct itself after some iterations. One solution to reduce the dependence on initialisation is to proceed by blocks; instead of updating π after each observation i, π could be updated for instance every 10 observations. Finally, we have to ensure that π sums to one, so in each step, it is divided by its sum. This is also useful, as it will never set any of its value to zero exactly. This means that any candidate variable still has a chance to be selected in the subset. Therefore, I propose a second version of the LLMS algorithm. The inclusion probability P depends strongly on the amount of penalty applied θ. Ideally, we would like to find the optimal amount of penalty θ that minimises the difference between S * and Ŝ.

From a simple look at figure 4, this difference is minimized when the discrepancy among the probabilities of inclusion is at its maximum. This result suggests the following optimal rule:

θ * = argmax θ (V ar [P θ ]) (33)
Figure 4 illustrates the selection paths P as well as the optimal rule for the penalty parameter. Also, when θ → 0 there is no penalty, and thus all the p candidate variables are selected:

P θ→0 = 1.
When θ → ∞ there is infinite penalty, and thus no candidate variables are selected: P θ→∞ = 0. This implies: V ar [P θ=0 ] = V ar [P θ=∞ ] = 0. However for intermediate amount of penalty 0 < θ < ∞ we have V ar [P θ=∞ ] ≥ 0. As P is decreasing in θ there is a unique maximum.

2.4.2. Reducing the sensitivity to the neighbourhood's size parameter via the q-kNN Algorithm

It is known that kNN is sensitive to the hyper-parameter k [START_REF] Sun | An adaptive k-nearest neighbor algorithm[END_REF].

Intuitively, the optimal value of k would be the one that produces the largest neighbourhoods, such that the linear assumption holds. It follows that this optimal value is dependent on the shape of the manifold, which is unknown. But even in the case I could measure distances on the R m space, these are still not the same as the ones on M. Indeed, take for example Fig. 5. Illustration of the revisited q-kNN algorithm.

(a) kNN (b) q-kNN Notes: The aim of both algorithms is to find neighbourhood (red dots) around a given data point (green dot): kNN (left) is using simple euclidean distances while the revisited q-kNN (right) is using compounded euclidean distances. The neighbourhoods both contain k points in the end.

The q-kNN comes at the cost of computations. For each neighbourhood, distances w.r.t to all the points have to be evaluated, several times. This can be shortened by computing all pairwise distances once. Then the remaining computations are only orderings, which are usually quite fast. However, if q-kNN is used at the same with time RkNN, it is worse. The pairwise distances have to be evaluated at each step, thus it is very slow.

Selection Properties

In this part, I study some properties of the selection estimator σ LLM S (X, f •g). The usual selection consistency property is known in the literature as the "Oracle Property" [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF]. The weak version of the Oracle Property is stated as follow:

lim n→∞ P ( Ŝ = S * ) = 1. ( 35
)
The strong version requires also the convergence of the estimated model, as established by [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF]. In that case, that would mean convergence of f • ĝ. But these are never explicitly estimated. Also, since I am primarily interested in the selected set, I will only focus on the Weak Oracle Property. Nevertheless, the Oracle Property in this form is not easy to work with. So first, I define the boolean version of the set of active variables S * denoted Ψ * . It is a vector of with p elements.

S * := {j : ∃j = j, ∂X j ∂X j = 0} ∈ N Ψ * = 1 j∈S * ∈ {0, 1} p . (36)
Their empirical local counterparts are defined as:

ŜN i := {j : ∃j = j, ∂X j ∂X j = 0, X ∈ N i } ∈ N ΨN i = 1 j∈ ŜN i ∈ {0, 1} p . ( 37
)
The Local Linear Manifold Selection operator, denoted P, is simply the average of the local selection operators after the boolean transformation:

P = 1 n n i=1 ΨN i ∈ [0, 1] p .
Now, the Weak Oracle Property can be restated, so it is easier to work with the LLMS estimator. lim

n→∞ P ( Ŝ = S * ) = 1 ⇔ lim n→∞ P ( Ψ -Ψ * 2 2 = 0) = 1 ⇔ lim n→∞ P ( Ψ -Ψ * 2 2 ≤ ν) = 1 ∀ν ≥ 0. (38)
I will refer to 38 as the Continuous Weak Oracle Property. Then, the LLMS estimator is consistent in selection if:

lim n→∞ P ( P -Ψ * 2 2 ≤ ν) = 1 ∀ν ≥ 0. ( 39
)
However, since P is the average of a linear selection operator Ψ N i , its selection properties have to be inherited from the selection properties of Ψ N i . The properties of the latter may vary substantially from one estimator to the other. So, I cannot directly deduce the properties of P in a general way. Nevertheless, the continuous weak oracle is illustrated using simulations in section 4.3. The results seem to indicate that the property holds.

Still, there is one thing of interest, which is to show the behaviour of P w.r.t to any given local estimator. Therefore, I make three propositions to show that P is never worse, in terms of selection consistency, than its local linear estimator Ψ N i . Following these results, by choosing a consistent local selection estimator ensures the selection consistency of the LLMS estimator.

Proposition 2. The LLMS estimator P is closer to the true set Ψ * than the underlying local linear selection estimator Ψ N i .

P ( P -Ψ * 2 2 ≤ ν) ≥ P 1 n n i=1 ΨN i -Ψ * 2 2 ≤ ν (40)
Proposition 3. When the level of type I and type II error α tends to zero, the underlying local linear selection estimator Ψ N i is closer to the true set Ψ * than the LLMS estimator P.

When α = 0, both are perfectly consistent.

lim α→0 P ΨN i -Ψ * 2 2 ≤ P -Ψ * 2 2 → 1 (41)
Proposition 4. When the number of variables increases, the LLMS estimator P is closer to the true set Ψ * than the underlying local linear estimator Ψ N i .

lim p→∞ P ΨN i -Ψ * 2 2 ≤ P -Ψ * 2 2 → 0 (42)
Proof of Proposition 2 is pretty simple, following the triangle inequality. First, I expand each of the norms.

P -Ψ * 2 2 = 1 n n i=1 ΨN i -Ψ * 2 2 1 n n i=1 ΨN i -Ψ * 2 2 = 1 n n i=1 ΨN i -Ψ * 2 = 1 n n i=1 ΨN i 2 + (Ψ * ) 2 -2 1 n n i=1 ΨN i Ψ * 1 n n i=1 ΨN i -Ψ * 2 2 = 1 n n i=1 ΨN i -Ψ * 2 = 1 n n i=1 Ψ2 N i + n(Ψ * ) 2 - 1 n n i=1 2 ΨN i Ψ * (43)
Then, I compare each terms of the expansions. The first terms are easy to compare by Jensen inequality

1 n n i=1 Ψ2 N i ≥ 1 n n i=1 ΨN i 2 .
The term in the middle is bigger if the sample size is greater than one.

n(Ψ * ) 2 > (Ψ * ) 2 ∀n > 1.
And the last terms in 43 are equal, so that concludes the proof.

Proposition 3 shows the behaviour of the LLMS estimator w.r.t the "consistency" of the underlying linear selection operator. That consistency is represented by a parameter α ∈ [0, 1], which can be understood as a Type I/II error. To do so, I need to assume that each element of the boolean linear selection operator ΨN i is a Bernoulli random variable.

π ≡ (π 1 , ..., π p ) ∈ [0, 1] p , π j =    1 -α if j ∈ S * α if j / ∈ S * ΨN i ∼ Bern(π). (44)
Since the vector Ψ * is only composed of ones and zeros, then it follows that the difference ( ΨN i -Ψ * ) is also a vector of Bernoulli r.v. having all the same parameter equal to α.

( ΨN i -Ψ * ) ∼ Bern(α). (45)
Also, because the linear selection vector is a boolean vector, taking its square does not change anything, it remains a boolean vector. Therefore its distribution is unchanged.

ΨN i ∈ {0, 1} p ⇒ ( ΨN i -Ψ * ) 2 ∼ Bern(α). ( 46
)
The first norm ΨN i -Ψ * 2 2 we are interested in, is the sum of the squared difference vector above. Because each element of the difference vector is a Bernoulli r.v., its sum, denoted y, is therefore a Binomial r.v. with parameter p, which is the length of the sum.

y ≡ 1 p p j=1 ( ΨN i -Ψ * ) 2 ∼ Binom(p, α).
Now, moving on to the second norm P -Ψ * 2 2 . By definition, the expected value of the local linear selection is the vector of parameters π in equation 44. It happens to be exactly the definition of the LLMS estimator.

E[ ΨN

i ] = π lim n→∞ P → π
We can exploit this result to obtain the norm in terms of α and p only:

lim n→∞ P -Ψ * 2 2 = π -Ψ * 2 2 = α 2 2 = p j=1 α 2 = α 2 p.
Knowing this, I can write the probability of the first norm, represented by the r.v. y, being greater than the second norm which is equal to α 2 p. That probability is easy to obtain as y is binomial.

P y ≤ α 2 p = α 2 p l=0 p l α l (1 -α) p-l
This probability tends to one as α decreases, and that concludes the proof. However it is worth noticing that when α = 0 then both norms are equal to zero.

Proof of Proposition 4, follows directly from the previous proof. Simply by taking the limit of the same probability as the number of candidate variables p increases.

lim p→∞ P Y ≤ α 2 p → 0.
That limits tends to zero, which means that the local linear selection is always worse than LLMS when the number of variable is very large. The LLMS smooths by averaging, and so remains closer to the true set of selected variables.

Global Approach via Diagonal Auto-Encoders

Compared to the first section, here I use a completely different approach. Instead of a local estimator, Neural Networks provide a global estimator of the manifold.

There have been plenty of papers on neural networks that implemented variable selection.

Many empirical works only uses some heuristics, while the theoretical part is significantly smaller. This literature is mainly driven by computer scientist, so here I also want to highlight the connection with more traditional statistics. Frequently, in regression models for instance, variable/feature is performed through constrained optimization. Parameters associated with a given variable are penalized such that some of them are exactly zeros, thus performing variable selection. These are known as "penalty methods", see Fan and Lv. ( 2010 Parameter penalization is quite easy to solve in the linear models. Whereas there is often no analytical formulas, simple algorithms can be implemented to get a solution. In non-linear models, the usual practice is to rewrite the problem using basis transformations [START_REF] Wang | Shrinkage estimation of the varying coefficient model[END_REF][START_REF] Ni | Automatic model selection for partially linear models[END_REF] to reduce it back to a linear problem, solved with almost the same aforementioned algorithms. A Neural network can indeed be subject to these kinds of constraint since a network is a combination of multiple basis transformations.

These bases are formed with non-linearities, often a sigmoid/tanh function, and parameters matrices. But only imposing sparsity on these matrices does not ensure variable selection, because a given variable is connected to multiple bases. So a sparse neural network is not necessarily a selection operator. One has to impose a sparsity condition on every connection of a variable to the network. Again, there exist different ways to implement this solution.

But before going into the final manifold selection operator, I will introduce a simple concept for sparse neural network, that ensures variable selection in the end.

Restricted Diagonal Layer for Variable Selection

First, I introduce the mechanism I use for variable selection in the neural network framework. To do so, I consider the basic non-parametric model, that distinguishes between inputs and outputs. The objective is still to find the sets of "true" variables, but now there are two such sets: S * x and S * y .

y = m(x) (X 1 , X 2 , ..., X px ) ≡ X ∈ R n×px (Y 1 , Y 2 , ..., Y py ) ≡ Y ∈ R n×py m : R px → R py S * x := {j : ∃j , ∂Y j ∂X j = 0} S * y := {j : ∃j , ∂Y j ∂X j = 0} |S * x | ≤ p x |S * y | ≤ p y ( 47
)
The objective is to build a sparse estimator that only include the true sets of variables.

And the mapping m is approximated by a neural network. For the choice of the structure of the neural network, I will stick to the proof of [START_REF] Hornik | Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks[END_REF] and use a 3-layer architecture (as illustrated in Figure 6). To summarise, it is the minimal architecture required to approximate any non-linear function. From the viewpoint of these authors, there is no need to make the network deeper, i.e. having more layers.

Fig. 6. Basic architecture for a Neural Network

X 1 X 2 X 3 Z (1) 1 Z (1) 2 Z (1) 3 Z (1) 4 Z (1) 5 Y 1 Y 2 W 1 W 2
For the sake of simplicity, I omit biases in the representation even though they are included at every layer, but said otherwise. Nevertheless, the simple 3-layer network does not provide a sparse solution. So how do we implement variable selection in neural networks?

In the Machine Learning literature, selection often consists in quantifying the predictive importance of a given variable/feature [START_REF] Breiman | Random forests[END_REF]). This can be very computationally consuming. To make the network select variables, I use a different approach based on [START_REF] Neuneier | How to train neural networks[END_REF] and [START_REF] Vecoven | Feature selection with deep neural networks[END_REF]. The main idea is to introduce diagonal matrices that directly connect inputs and outputs to the network, controlling their entry in the model. W.r.t these papers, there are two contributions: (1) selection is performed also on the outputs, not only the inputs, (2) the domain of the diagonal selection parameters is constrained to be between 0 and 1. Figure 7, depicts the addition of the two selection layers at both ends of the network. On top of this, I do not need to enlarge the loss function to impose restrictions on parameters, contrary to [START_REF] Ye | Variable selection via penalized neural network: a drop-out-one loss approach[END_REF] for instance. Only the structure of the network is modified. Restricting the selection parameters to the [0, 1] domain has advantages in terms of interpretation: they can be seen as "pseudo-probabilities".

If a variable seems to be useful to the network then its associated parameter in the input selection layer tends to 1, otherwise, it tends to 0. That directly links to the LLMS estimator presented in the previous part of the paper. Both give a continuous selection vector, whose elements are between 0 and 1.

Fig. 7. Modified architecture for selection in Neural Network X 1 X 2 X 3 X1 X2 X3 Z (1) 1 Z (1) 2 Z
(1) 3

Z (1) 4 Z (1) 5 Ỹ1 Ỹ2 Y 1 Y 2 W s 1 W 1 W 2 W s 2
I restrict the domain of W s 1 and W s 2 using the non-linear sigmoid function.

sigmoid : R → [0, 1] γ → 1 1 + e -γ . W s 1 = diag [sigmoid(γ 1 )] , W s 2 = diag [sigmoid(γ 2 )] .
For the rest of the network, I use the hyperbolic tangent as the activation function, mainly for the reason that it is symmetric. The output layer is linear so that the network can map any real inputs to any real outputs, without further restrictions. For optimization stability one may still normalize the variables w.r.t their means and standard deviations.

No activation function is applied to the selection layers. The full set of equations, in matrix form, for the modified network structure is:

W s 1 = diag [sigmoid(γ 1 )] X = XW s 1 Z (0) ≡ X Z (l) = tanh(Z (l-1) W l ), ∀l ∈ 1, • • • , n h -1 Ŷ = Z (n h -1) W n h Ỹ = ŶW s 2 W s 2 = diag [sigmoid(γ 2 )] E = Y - Ỹ
Where n h represents the depth of the network, the number of successive hidden layers. It is a hyper-parameter, that is chosen arbitrarily. Its value is in fact very linked to the neighbourhood size k in LLMS. Both control the degree of non-linearity, and ultimately depends on the shape of the DGP.

The optimization program is: min

W 1 ,...,Wn h ,γ 1 ,γ 2 T race(E E).
The estimated selections are represented by selection weights W s 1 and W s 2 obtained at the end of the optimization:

Ŝx → diag(W s 1 ), Ŝy → diag(W s 2 ).
To illustrate the validity of the concept, I run a simulation according to model 47. The selection paths are reported in Figure 8. Green paths correspond to the variables that belong to the true sets S * x and S * y , the irrelevant ones are coloured in red. The green paths converge to 1 and the red ones to 0, so it is quite obvious from the chart which variables should be selected. Therefore, the selection seems consistent is this simulation. Nonetheless, I observe the paths for inputs and the paths for outputs do not converge in the same way. On the one hand, the inputs' parameters spread from epoch 1 to epoch 2000, and then are attracted to the boundaries. On the other hand, the outputs' parameters all go down first, and then diverge around epoch 2000.

Diagonal Auto-encoders Manifold Selection

In this section, I extend the model from equation 47 to take into accounts non-linear manifolds. It is a generalisation, and it encompasses the previous section as a special case.

Inputs X and outputs Y are embedded into a the larger "new inputs" X of size p = p x + p y .

So the network has the same data at both ends. The model is therefore the same as for LLMS:

M :={x ∈ R p |x = f • g(x)} g : R p → R r , f : R r → R p , p > r (X 1 , X 2 , ..., X p ) ≡ X ∈ R n×p S * := {j : ∃j = j, ∂X j ∂X j = 0} |S * | ≤ p (48)
The projection function g maps the input onto a low-dimensional space, that is mapped back to the original space by the immersion function f . The problem turns into the estimation of two separate non-linear functions. Each function can be approximated using the 3-layers neural network. So to estimate both functions, the network is enlarged by combining two 3-layers architectures. At the center, they are connected by a "contraction" layer. This layer reduces the dimension of the input to the dimension of the manifold, as depicted in Figure 9. Its dimension has to be lower than the one of the input, otherwise, the network may learn the identity function. This architecture is not new. It is known as the 5-layers Auto-Encoder, proposed by [START_REF] Kramer | Nonlinear principal component analysis using autoassociative neural networks[END_REF].

As I did before, I add the selection layers to this Auto-Encoder, thus obtaining a 7-layers architecture. Since the input flow through a deeper network, I will refer to this selection estimator as DAMS, for Diagonal Auto-Encoder Manifold Selection. Also, because the inputs and outputs are now the same, I can enforce the two selection parameters to be the same.

Thanks to the chain rule, the derivative of this shared weight is fairly easy to compute. I only have to sum the two derivatives.

process. The nice thing is that I can still provide the so-called "selection paths" for DAMS exactly as for LLMS. It suffices to plot the evolution of the parameters diag(W s ) throughout the optimisation process. Obviously, the optimal value for W s is chosen to be the last of the optimisation process.

Smoothing selection paths with Ensemble Averaging

There are two reasons why the paths of the parameters might be unstable. First, optimisation of both selection weights and the rest of the network parameters simultaneously. This means that the network has to learn the manifold at the same time it does the selection.

A better idea would be to start learning the manifold with all the variables, even though it means a less precise estimator. After that, start the selection. The estimator will then be more precise. So W 1 , ..., W n h are being optimised few epochs before W s . Second, the parameters estimated for a given network depend on weights initialization. A single set of selection paths W s might not be stable, nor lead to consistent selection. This phenomenon is depicted in Figure 10. The green paths correspond to the variables that are in S * , while the red correspond to the noisy variables. Indeed, the paths are sometimes chaotic. The solution is to average the paths across many random initialisations. This method is known as "Ensemble Averaging". The network is replicated a certain number of times, to form what is called an "Ensemble". The paths for every replication of the network are then averaged.

That produces more stable, but also more consistent selection paths, as it is shown in Figure For each simulated process, A is randomly chosen, so that the shape of the manifold changes as well. Very small changes in A might induces very different manifolds in the end. I simulate

X j = f j (Z j A j ) with Z j ∼ U [-2, 2], j = j , A ∼ U [-2, 2].
The non-linear transformations are chosen among the following ones:

• f 1 (x) = x, • f 2 (x) = x 2 , • f 3 (x) = sin x 2 , • f 4 (x) = cos x 2 ,
• f 5 (x) = tanh (2x),

• f 6 (x) = exp (x-10) 2 20 , • f 7 (x) = (x + 1)1 x>x -(x -3)1 x≤x .
I also set the following parameters:

• n = 5000,

• p = 50,

• d = 7,

• υ = {0.01, 0.25}

There are two hyper-parameters left. The size of the neighbourhood k for LLMS, and the size of the hidden layers for DAMS. Both control the degree of non-linearities.

I choose to perform both algorithms with 3 different values of these parameters. For LLMS, k = {5%, 20%, 50%} of the sample size n. And for DAMS, the hidden layer sizes are {(7, 1, 7), (7, 3, 7), (7, 5, 7)}. The simulation is run 1000 times and results are averaged. The computations are performed on an Intel(R) Core(TM) i5-3220M CPU 2 × 2.60GHz with 12GB RAM. The estimated selections are compared with existing methods, both parametric and non-parametric. I only compare to selection methods which do not make difference between explanatory variables and outcome variables. These are methods which act on covariance/correlation matrices, for instance, the LASSO is not valid for comparison in these simulations but the Graphical LASSO is.

• Eigen Thresholding (ET)

Σ = VLV p (Σ, θ) = V ij 1 |V ij L 1/2 j |>θ
• Sparse PCA (sPCA) [START_REF] Zou | Sparse principal component analysis[END_REF] Σ = (XVV ) (XVV )

p (Σ, θ) = θ 1 V 1 + θ 2 V 2 2 • Graphical LASSO (gLASSO)
p (Σ, θ) = θ 1 U 1 + θ 2 V 1 • NonParanormal (NPN) Lafferty et al. (2012) f X = (f 1 (X 1 ), ..., f p (X p )) f j = E(X j ) + V ar(X j )Φ -1 Fj (X j ) ∀j ∈ {1, ..., p} Σ = Cov(f X ) p (Σ, θ) = log(det(Σ -1 )) -tr(Σ -1 ) -θ Σ -1 1
where Φ -1 is the Inverse Normal CDF and F is the empirical CDF. Results for low variance of noise are reported in Table 1. Each column corresponds to a different DGP. There are six possibles DGP, three linear, three non-linear. Each has a different reduced dimension of r. When that dimension is equal to 1, visual inspection is enough to assess the correlation, as you can see in the graph d) from Figure 12. The correlations become less obvious as r increases. Therefore, the value of r may have an impaction on the selection properties. So in the end, there are four columns in the tables. The first one aggregates the three linear models since the results I obtain do not vary at all with r. Then, each line of the table corresponds to a different algorithm. The three first blocks are for the three versions of LLMS. So there are 9 lines corresponding to LLMS, since there are three values for k. The next block is for DAMS, there are also 3 lines. The last block is for all other algorithms.

In the first column, you observe TPR (left) and FPR (right) for linear models. Without any surprises, all algorithms have 100% performance. They all select the true set S * and have also zero FPR.

Next, in the second column, there are results for the non-linear manifold of dimension r = 1.

Overall, LLMS performs very well. I observe that LLMS3 is very slightly superior to LLMS1.

This result shows that correcting distances with q -kN N improves selection. Another finding is the choice of k does not seem to greatly affect the TPR. That finding is also observed for DAMS, whatever the hidden layer sizes it always selects the true set. Surprisingly, the other algorithms also have decent performances on the non-linear DGP. The linear selection operators (ET, gLASSO and PMD) have a near 90% TPR. Only PMD has a non-zero FPR, but it is quite low (7.2%). The only other non-linear operator which is NPN, has 90.3% TPR and 5.7% FPR. I can conclude that in case of very low dimensional non-linear manifold (r = 1), the linear operators have quite good selection properties.

Finally, in the last columns (r = 2 and r = 3) the performances change radically. The sensitivity to k in LLMS is much greater, whereas DAMS is robust to different layer sizes. This result highlights the difference between a local and a global approach. If k is at 50%, all versions of LLMS fail to select the true set. Also, LLMS1 and LLMS3 have much poorer performance than LLMS2, respectively 85.7% and 86.3% versus 100%. Among the proposed methods, only DAMS with the largest hidden layer (7, 5, 7) has non-zero FPR. This is due to the fact that, in these cases, the contraction layer has a dimension greater than r. Therefore, there are too many dimensions than what is needed to approximate the manifold. These extra dimensions are used by the network and learn the identity function for some noise variables. Symmetrically, the smallest structure (7, 1, 7) has only 85.7% TPR. There are not enough dimensions to approximate the manifold, and it biases the selection afterwards. So, just as k, there is a trade-off for choosing a good value of the hidden layer sizes. Nevertheless, the performance is better than the other algorithms. The linear operators exhibit two patterns. The first one is characterized by a low TPR (around 20%) and a zero FPR. The second pattern is characterized by a relatively high FPR (around 70%) but with also a high FPR (ranging from 40% to 80%). So in the first case, they fail to select all the relevant variables, only some of them. In the second case, they simply select too many variables, including, by chance, some of the relevant ones.

The runtimes show that both LLMS and DAMS are much longer to perform than any competitors. Among the three versions of LLMS, LLMS3 is the longest. The best version in terms of TPR, is LLMS2. It takes around a minute, which is reasonable. The computations depend on the neighbourhood size k, so the runtime can be reduced by reducing k but that may affect selection. DAMS takes, on average, 6 minutes. The computations heavily depend on the batch size and the number of networks in the ensemble. Reducing these two parameters can shorten the runtime but may affect selection. 1 is the level of the noise. This noise may affect the selection, as it is more difficult to approximate the manifold.

In the first column, for the linear DGPs, I observe almost the same performance for all algorithms, above 95%. DAMS is the only algorithm to have 100% TPR. In columns 2, 3 and 4, the overall performance is poorer than in the low noise case. The relative performance each value of ν.

The results of the simulations for LLMS2 are presented in Figure 16. In the upper part, the norm decreases with sample size. In the lower part, the observed frequency tends to 1 as n increases for any value of ν. This finding supports the hypothesis that LLMS is indeed an Oracle estimator. I also note that the convergence is pretty fast. Most of the coloured lines are already closed to 1 at sample size 5 000.

The results for DAMS are presented in Figure 17. In the upper part, the norm decreases with sample size as well. However, for a small value of n, there are some outliers. The DAMS operator is therefore not stable in selection, in small samples, in these simulations.

As n increases, the observed frequency tends also to 1 for any value of ν. Again, this result supports the hypothesis that DAMS is indeed an Oracle estimator. Nevertheless, the convergence is slower than LLMS. Still, there is room for improvement, especially for hyper-parameter selection. For LLMS, I propose an optimal rule for the penalty parameter θ, but not for the neighbourhoods' size k. For instance, the latter could be optimized w.r.t to the manifold reconstruction error.

For DAMS, the only hyper-parameters of interest are the layer size, especially the central contraction layer. If it is larger than r, there can be over-identification of f and g. But it is not obvious that selection parameter will be affected. If the contraction layer size is larger or equal to |S * |, the network can start to learn the identity function, and select more variables.

Unfortunately, both r and |S *

x | are unknown beforehand. These are issues for future research to explore.

)Fig. 2 .

 2 Fig. 2. Local linear approximation of a non-linear manifold at t 0 = π 4 , -π 4 , 3π 4 , -3π 4

I

  perform non-parametric variable selection on manifolds using a linear selection operator locally, i.e. on neighbourhoods of the data. Solving the selection at the global level is hard because of the unknown possible non-linearities. However, simple linear selection operators have been well studied and greatly improved over time. The major restriction they suffer is the linearity assumption. To ensure linearity holds, I slice the global non-linear manifold into lots of neighbourhoods. Inside a given neighbourhood, a simple linear selection operator should work very well. Doing so produces a selection Ŝ for each neighbourhood. As I want a single selection in the end, I just merge all the selections by averaging them. Therefore, the original problem of performing selection on non-linear manifold is solved.

Fig. 3 .

 3 Fig. 3. Example of non-linear DGP and LLMS selection paths

Proposition 1 .

 1 If S * ∈ B and m ≡ |B| < |A| = p then RkNN yields closer neighbourhoods on M than kNN.

Algorithm 2

 2 Local Linear Manifold Selection with random subsets (LLMS2) input X, k, Θ, σ(X, Ω)Start with equal subset probability π = 1 p for i = 1 to n do Choose a subset of {1, ..., p} from probability the probability vector π Search for the neighbourhood N i,k (X i ), computing distances w.r.t to the subset.for θ ∈ Θ do Estimate the sparse correlation matrix RN i,k using σ(N i,k , Ω) Record active variables in ŜN i,k ,θ = j| RN i,k (j, s) = 0, ∀s = j . Turn Ŝ into boolean vectors Ψ N i,k ,θ = 1 j∈ ŜN i,k ,θend for Update the subset probability vector π = π + 1 i,k ,θ . output P θ 2.4. Optimal Choice for Hyper-parameters 2.4.1. Choice of penalty parameter

Fig. 4 .

 4 Fig. 4. Choice for optimal penalty

  ); Huang et al. (2012); Mehmood et al. (2012); Jović et al. (2015); Desboulets (2018) for complete reviews.

Fig. 8 .

 8 Fig. 8. Selection paths for inputs and outputs selection

3. 3 .

 3 Estimation of the Parameters of the Network 3.3.1. Stochastic First-Order Optimization The architecture proposed only enlarge the network by 2 linear layers. I do not modify the loss function. Thanks to the chain rule, it makes the computations of the derivatives for each weight as simple as any other layers. Standard backpropagation algorithm can be used without further modifications. Thus, I simply use the gradient descent algorithm, computed on a random batch of the data. This method is known as stochastic gradient descent. It has two main advantages: (1) it avoids local minimas more easily and (2) it reduces the duration of the computations. The size of the batches is a hyper-parameter, usually set to a fraction of the total number of observations.

Fig. 10 .

 10 Fig. 10. Examples of DAMS selection paths with a single network

1 •

 1 [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] Σ = Cov(X) p (Σ, θ) = log(det(Σ -1 )) -tr(Σ -1 ) -θ Σ -1 Penalized Matrix Decomposition (PMD)[START_REF] Witten | A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis[END_REF] Σ = (dUV ) (dUV )

Fig. 12 .

 12 Fig. 12. Simulated DGP with low noise (υ = 0.01)

Fig. 16 .

 16 Fig. 16. Continuous Oracle property -LLMS

  Third, this implies that all candidates are measured with proportional noise V ar [ j ] = υV ar [X j ]. If I were to use a covariance matrix, this would imply constant variance of the noise: V ar [ j ] = υ. However, keep in mind that this assumption holds only locally. It can vary across the sample, as we compute it in different neighbourhoods.

Table 1 :

 1 Selection rates -Simulation with low noise (υ = 0.01)

		Linear		Non-Linear
		r = {1, 2, 3}	r = 1	r = 2	r = 3
		TPR FPR	TPR FPR TPR FPR TPR FPR Runtime (sec.)
	LLMS1 k=5% 100. -LLMS1 k=20% 100. -LLMS1 k=50% 100. -	99.9 -98.9 -96.4 -	88.9 -83.5 -59.4 -	85.7 -71.4 -51.9 -	61.200 97.300 291.800
	LLMS2 k=5% 100. -LLMS2 k=20% 100. -LLMS2 k=50% 100. -	100. -100. -96.5 -	100. -98.6 -70.6 -	100. -98.7 -71.1 -	50.800 60.800 73.400
	LLMS3 k=5% 100. -LLMS3 k=20% 100. -LLMS3 k=50% 100. -	100. -100. -96.8 -	90.6 -85.2 -62.1 -	86.3 -73.4 -50.3 -	179.300 306.500 1775.300
	DAMS (7,1,7) 100. -	100. -	100. -	85.7 -	342.700
	DAMS (7,3,7) 100. -	100. -	100. -	100. -0.1	354.800
	DAMS (7,5,7) 100. -	100. -	100. -4.6 100. -9.8	391.200
	ET	100. -	86.6 -	37.1 -	17.6 -	0.014
	sPCA	100. -	23.2 -	30.3 -	25.8 -	0.625
	gLASSO	100. -	89.3 -	46.7 -	18.7 -	0.036
	PMD	100. -1.9	89.6 -7.2	87.5 -25.6	70.8 -40.1	2.279
	NPN	100. -	90.3 -5.7	70.9 -28.9	71.4 -76.6	1.058
	All values equal to zero are left as blanks.

Table 2 :

 2 Selection rates -Simulation with high noise (υ = 0.25) All values equal to zero are left as blanks.Results for high variance of noise are reported in Table2. The structure of the table is the same as the preceding. The only difference w.r.t to Table

		Linear		Non-Linear
		r = {1, 2, 3}	r = 1	r = 2	r = 3
		TPR FPR	TPR FPR TPR FPR TPR FPR Runtime (sec.)
	99.1 -LLMS1 k=20% 99.3 -LLMS1 k=5% LLMS1 k=50% 97.7 -	94.3 -88. -83.4 -	90.4 -73.7 -55.8 -	76.4 -57.1 -41.9 -	59.600 104.100 313.300
	99.7 -LLMS2 k=20% 99.7 -LLMS2 k=5% LLMS2 k=50% 99.1 -	99.7 -97.3 -85.4 -	99.8 -98.1 -68.5 -	99.6 -95.7 -50.6 -	51.000 60.000 69.100
	99.4 -LLMS3 k=20% 99.3 -LLMS3 k=5% LLMS3 k=50% 98.0 -	97.1 -90.9 -84.9 -	100. -75.6 -59.4 -	80.7 -58.4 -41.2 -	147.400 300.400 1825.300
	DAMS (7,1,7)	95.7 -	94.4 -	92.3 -	74.7 -	334.700
	DAMS (7,3,7) 100. -	100. -	100. -	97.1 -	355.300
	DAMS (7,5,7) 100. -	100. -	100. -	100. -	371.200
	ET	93.7 -	76.6 -	44.3 -	27.5 -	0.016
	sPCA	97.6 -	40.6 -	32.3 -	27.1 -	0.635
	gLASSO	97.4 -	80.3 -	45.1 -	24.6 -	0.030
	PMD	96.7 -5.8	86.9 -18.6	85.9 -28.5	82.9 -35.5	2.170
	NPN	99.9 -0.1	82.6 -	71.4 -35.7	72.9 -73.2	1.093

the very usual Euclidean distance. If the manifold is very non-linear, increasing euclidean distances will produce neighbourhoods in which points are very apart on the manifold (see Figure 5(b)). However, small euclidean distances approximate well the distances on the manifold. I can use this property to construct better neighbourhoods. This does not solve the problem of the optimal value of k but still reduces the sensitivity of the algorithm in case the manifold is very non-linear.

So I propose a refined version of kNN, in which neighbourhoods are constructed sequentially from small distances only. The basic idea is to take the q < k nearest points to x i and then expanding along with the next q nearest points to the ones previously selected, etc... until the neighbourhood's size is k. This better ensures the local linearity assumption and tends to reduce the sensitivity to k. Therefore I call the new version q-kNN.

The first neighbourhood contains q points. Then I compute the q-nearest points of these q points, so the second neighbourhoods as q 2 points. However doubles are removed, so I end up with a second neighbourhood having between q + 1 and q 2 points. So in the n th step the n th neighbourhood of size at least q + n and at most q n . With a not too small value of q this sequence quickly converges to a neighbourhood with k points.

The full set of equations, in matrix form, for the DAMS structure is almost the same.

The optimization program is:

And the estimated set is simply:

The DAMS estimator diag(W s ) is an analogue to the LLMS estimator P. Their main difference lies in the differentiability and continuity assumptions of the manifold. DAMS requires both f and g to be continuous everywhere, while LLMS requires continuity only locally, inside the neighbourhoods. Also, for LLMS, the two functions have to be differentiable, at least locally. Differentiability is not needed for DAMS. If the true manifold is not continuous everywhere, DAMS might be inconsistent. The selection parameters depend on the derivatives of the total error. They are not linked to the derivatives of f or g, contrary to LLMS.

Another thing is that the LLMS estimator is dependent on a penalization parameter θ.

That is not the case for DAMS, which incorporates that mechanism inside the optimisation between each algorithm and the others stays the same, compared to the previous table. Only LLMS2 and DAMS are consistent estimators for all values of r. Still, LLMS2 is sensitive to k. A value between 5% and 20% produces high TPR, ranging from 95.7% to 99.8%. For DAMS, the (7, 5, 7) structure always produces perfect selection.

The results from these simulations are the following. LLMS2 and DAMS are good candidates as selection estimator for data on non-linear manifolds. The basic implementation (LLMS1) is not consistent enough, and it is very sensitive to k. The last version LLMS3 is indeed better than LLMS1, but not as much as LLMS2 is. A combination of LLMS2 and LLMS3 (r-kNN + q-kNN) would certainly produce even better results. But it would need even more computations. And when looking at runtimes of LLMS3, such a combination would be computationally prohibitive.

A linear selection estimator fails to select the right set of variables when the true DGP is non-linear.

Numerical Oracle Properties

The Oracle Property is the property for a selection operators to uncover the set of true variables S * . This property has already been discussed for LLMS in section 2.5. Because the two proposed estimators take on continuous values, I have defined that property in terms of a distance w.r.t to the true set, in boolean form, denoted Ψ * . Therefore, the Weak Continuous

Oracle property for LLMS and DAMS is stated as:

with ν arbitrarily small. So in this section, I investigate on that property numerically for the two algorithms. For increasing values of n, I simulate a non-linear DGP 1000 times, as I did in the previous section. Then, I compute the frequency at which the observed norm is lower than an arbitrarily small value ν. If the estimator is consistent, as the sample size increases I should observe the norm decreasing as well. In the upper part of the figure, I report the norm along with sample size. It takes on four values n = {500, 5 000, 50 000, 500 000}.

The solid black line represents the median of the different runs for a given sample size. The lower part of the figure represents the frequency at which the norm is lower than a given small value ν. These small values range from 0.01 to 0.0001. There is a coloured line for


Natacha Aveline-Dubach

Covid-19 in China: the pandemic exacerbates the speculative mechanism in residential real estate

One of the striking effects of the Covid 19 pandemic has been the dramatic increase in the disconnection between finance and the real economy in the United States and Europe. In China there is a similarly growing disjuncture between house prices and economic fundamentals. Both cases show a process of asset inflation driven by a frenzy of investment in a context of deflationary pressure. But the comparison ends there, because the dynamics of real estate investment in China respond to the logic of a sui generis 1 capitalist regime that has placed the extraction of land value at the heart of its economic model.

), which local governments sought to close by commodifying public land through the use of their power to sell land use rights 2 . During the period of intense manufacturing production that made China the "world's factory", low-cost land use rights were granted in the urban periphery to manufacturing firms (initially foreign, then also domestic) with the aim of increasing employment and fiscal resources. Then, as industrial growth stimulated urbanisation, which generated higher land values, local governments sold land use rights in central and peri-urban areas at high prices to property developers (with little compensation from residents or farmers) for the construction of office and housing projects 3 . Urban infrastructures were financed by dedicated investment platforms (the famous LGFVs, Local Government Financing Vehicles) guaranteed by municipal land reserves, beyond the control of the financial authorities. These extra-budgetary revenues made it possible to cover the cost of urban growth while providing metropolitan areas with quality infrastructure, thus increasing their attractiveness for industrial production. In return, these strategies ensured an upward career for municipal executives (primarily mayors and party secretaries), under a system of executive promotion based on the ability to generate economic growth.

Figure 1. The growing budget deficit of local governments in China

Source : Natacha Aveline-Dubach, based on the China Statistic Yearbooks

Financing urbanization through households' real estate investments

However, for this virtuous circle to work, it was imperative that property developers successfully market their real estate operations downstream of the financing circuit. It was then mainly households that fed the 'land loop' (figure 2). Everything possible was done to encourage them to invest massively in the residential markets: dissemination of the property ideology as the basis of values centered on the family, development of financial circuits dedicated to home ownership, strong restrictions on financial investment. Several developers even started to sell commercial or hotel space to households on the basis of guaranteed income 4 . The adoption of an ambitious 4 trillion-yuan recovery plan after the 2008 crisis only accelerated this process by supplying conventional or informal real estate credit channels via investment platforms. Households therefore massively invested their savings in real estate, starting with housing markets whose extraordinary performance contributed to the selfvalidation of their speculative anticipations.

Then came the time when this situation became very cumbersome for the Chinese authorities. In the mid-2010's, housing prices became totally unaffordable for first-time buyers in most major cities, while the suburbs were surrounded by "ghost" neighbourhoods of vacant highrise buildings owned by multiple-owner households. The increase in housing oversupply has reached such proportions that the state has dispensed with any official assessment. The latest estimate, which dates back to 2013 and is very likely to be underestimated, shows 52 million empty flats, or 22% of the residential area in urban China (Figure 3). While most of this stock In the United States or in Europe, overproduction on such a scale, combined with a proliferation of private debt, would certainly have led to a fall in housing prices. In fact, the Western media alerted public opinion to the imminent deflation of the Chinese real estate bubble in the early 2010. However, quite the opposite happened during the following decade. Not only did house prices not fall, but they continued to rise, with an acceleration in the largest metropolises such as Shenzhen, Beijing and Shanghai (Figure 4).

Commercial real estate hit by the pandemic

While the lockdown froze flat sales for a few months, transactions gradually resumed without any negative effect on prices. However, the health crisis has had a serious impact on commercial property, where overcapacity was also evident at the end of 2019. According to converging estimates from the real estate industry, one fifth of A-grade office space is currently unoccupied in Shanghai and Shenzhen, and one third in deca-million cities such as Tianjin and Chongqing. There is no information on the occupancy rate of the less prestigious office buildings in suburban districts, but it can be assumed that it has deteriorated sharply. The effects of the health crisis have been, as elsewhere, even more severe in the retail and hotel sectors. The markets have taken note of these evolutions by lowering rents, and the real estate industry is organising itself to manage the surplus surface area by adapting it to postcovid practices: conversion of offices into spaces for flexible use, transformation of

2 0 0 1 2 0 0 2 2 0 0 3 2 0 0 4 2 0 0 5 2 0 0 6 2 0 0 7 2 0 0 8 2 0 0 9 2 0 1 0 2 0 1 1 2 0 1 2 2 0 1 3 2 0 1 4 2 0 1 5 2 0 1 6 2 0 1 7 2 0 1 8 Beijing Tianjin Shanghai Chongqing Shenzhen Guangzhou 35 villes CNY/m 2 7000 €/m 2 4730 €/m 2 3660 €/m 2 1580 €/m 2
cities commercial surfaces into logistics spaces and repositioning of tourist real estate towards services for senior citizens. The central government has made it clear that it will not come to the rescue of the real estate industry, but the local authorities are assisting developers with policies to support the so-called "first shop" economy (shoudianjin) aimed at replenishing shopping malls with new brands unknown in China 6 . Foreign investment, which had been virtually eliminated from the real estate market since 2008, is even being used to organise the debt reduction of Chinese companies through the sales of property assets 7 .

State efforts to curb the rise in residential prices

In the face of the turbulence that is shaking commercial real estate, the resilience of residential markets comes as a surprise. In Shenzhen and Beijing, prices have surpassed those in Hong Kong, averaging between 12,000 and 14,000 euros per square meter in the central districts, more than 40 times the local annual income 8 . The media are returning to the speculative bubble hypothesis, which has been put on the back burner in recent years. However, the central government has not spared its efforts to contain real estate inflation and absorb residential vacancies. It has tried to direct the flow of internal migrants towards smaller cities by delivering massively urban hukou -without much success -while real estate acquisitions in the largest cities have been limited to the holders of the local hukou 9 . Credit conditions have also been tightened: in hot markets, the amount of the initial deposit required per family has been increased to 50-70% for the purchase of a second home and to 100% for a third (in cities where such transactions are allowed), compared to 30% for a first purchase. Households quickly found the solution: divorce to benefit a second time from the conditions of the first purchase or raise funds on crowdfunding platforms to reach the required down payment amount 10 .

Private developers have been also targeted by state regulatory measures. They have been asked to redeploy their activity outside the Chinese residential market, at the risk otherwise of seeing their access to land dry up. The major property developers had already begun to expand their residential production activities abroad, offering their Chinese clients new investment opportunities in the United States, Canada and Australia, whose destabilising impacts on local markets had been in the news. But the offshoring of Chinese residential production is now taking a more spectacular turn in Southeast Asia, with the production of vast closed enclaves for Chinese residents, real pieces of largely vacant cities along the "Silk Road" corridor (Belt and Road Initiative project) 11 . In addition, awareness of China's ageing population has also enabled developers to obtain land to produce large residential developments for independent senior citizens and then market these homes to investors of all ages.

One way to curb the expansion of vacancies, or even reduce some of them, would be to develop private rental housing. This sector is highly devalued and poorly regulated in a homeowner society (about 84% of the urban population) where home ownership, the cornerstone of social welfare and a symbol of wealth, is affirmed even in marital practices 12 .

The State has made it a priority to build up quality private rental housing in the big cities.

Several large private developers have engaged in this sector with the aim of obtaining land for their other projects, but this activity is not compatible with the short-term nature of their financial resources -bank credit being restricted -and proves to be very unprofitable, with only 6% margin compared to 25% for housing sales in 2019 13 . It is therefore unlikely that this sector will see significant development in the current conditions of residential production.

The contradictory objectives of local policies

It must be said that local governments do not show a great appetite for change in their housing policy. Caught between the State's injunctions to curb the speculative mechanism and the need to finance productivist policies to ensure career advancement, they are struggling to break out of their addiction to land income. However, piecemeal rent capture through real estate projects is becoming more difficult to implement because of major changes at two major levels of the land loop : on the one hand, the increasing competence of the expropriated populations in the negotiation of compensation has considerably increased the cost of land resources, and on the other hand, the transformation of the productive system (development of highly innovative technologies and services, robotisation of factories and relocation of the manufacturing apparatus to South-East Asia) weakens the dynamics of urban growth, especially in the already overpopulated megacities.

In order to finance their economic and social policies, local governments have been allowed since 2015 to borrow on the bond markets. However, government-imposed bond quotas do not allow them to bridge the huge gap between their financing needs and their sources of income. As a result, they continue to use informal investment vehicles (LGFVs) for their infrastructure projects, even though the central government wants to put an end to this dangerous proliferation of "grey finance", which was estimated at between US$4,500 and 6,000 billion in 2018 14 . Municipal budgets could be supplemented by a more regular drain on wealth through property taxes, but local governments fear that the introduction of such a fiscal tool -a constantly-recurring issue-could lead to a fall in housing prices, and consequently a drop in the amount of land revenues 15 . The widespread adoption of a property tax would also be very unpopular in a society of homeowners, many of whom are heavily indebted and would not welcome the inequalities in wealth that a land registry could reveal.

As nothing can curb the inflationary dynamics of residential real estate, the State has called on local governments to develop an affordable homeownership offer. In the absence of any state financial support, the implementation of such housing projects is entrusted to property developers (private or public) on land that local governments are forced to concede at low cost. Unsurprisingly, this residential production remains very insufficient in view of the huge needs for affordable housing. Moreover, it only targets the holders of the local hukou, leaving aside the large mass of internal migrants who make up the bulk of the demand for low-cost housing. As a result, informal ownership rights known as "small property" (xiaochanquan) have developed, allowing "floating populations" to accumulate capital in urban real estate. In order to buy social peace, the authorities turn a blind eye to these spontaneous initiatives which have so far concerned some 250 million internal migrants (71 million households) 16 . No one knows how these rights will evolve, but the scale of the phenomenon is such that the State will probably not be able to continue to ignore them for long.

Increased polarization of real estate investments in major city-regions due to the pandemic

It is in this already largely problematic context that the Covid 19 pandemic occurred. The drop in sales during the lockdown period was only a brief episode in the residential markets, but it led to profound changes in the real estate industry. Small private developers, generally excluded from bank credit, did not survive the shock of the sales freeze. A concentration of the industry has therefore begun, with nationwide private developers managing to refinance their debt on the bond markets thanks to lower interest rates, and public developers (State-Owned Enterprises) being much better protected from the risk of bankruptcy due to their privileged links with the major public banks17 .

On the household side, the demand for investment in residential property in the largest cityregions is continuing, and even exacerbating. Investment flows tend to focus on first and second rank cities due to lower investment risk in vibrant economies and large residential markets. As in other countries, household investment in residential markets is driven by uncertainties about economic and financial developments. In the case of China, it is exacerbated by the need to build up precautionary savings to meet, inter alia, the costs of health, education, pensions and elderly care in a rapidly ageing society with residual social protection. This is all the more so as alternative investment vehicles are very poorly developed.

But what is most peculiar to China is the near-absolute confidence that homeowner households have in the government's ability to effectively steer residential markets [START_REF] Aveline-Dubach | China's housing booms: a challenge to bubble theory[END_REF] . There are several reasons for this : firstly, the very recent experience of a residential market that is only thirty years old and has not undergone any major turnaround to date 19 ; secondly, the close attention that the government pays to the evolution of property prices, knowing that a significant devaluation of this enormous savings reserve would be likely to seriously undermine the legitimacy of the current political regime; finally, the remarkable efficiency of the agile regulation of the markets by President Xi Jinping, a leader whose stability of power is assured for a long period of time. Such agility is permitted by the authoritarian nature of a regime in which the entire administration system is under the control of the Party. The state at its various levels can intervene "just in time" in the entire real estate value chain through its control of land resources and the banking system. The refusal to connect real estate and popular savings to global financial circuits, as demanded everywhere by Wall Street, has so far preserved the effectiveness of the Chinese party-state's levers of action, perhaps in return protecting the world economy from a major systemic crisis such as the GFC 20 . This is not without significant costs and risks: a huge waste of material resources embedded in deteriorating built-up areas -and contributing to the degradation of urban landscapes -an acute shortage of affordable housing in major cities, the proliferation of informal debts and property rights, and above all the risk that investors' shared belief in the government's regulatory capacity will be undermined. To alleviate the pressure on precautionary savings, it would be necessary to develop a social protection system that covers the needs of the population and puts an end to the blatant inequality between rural and urban hukou holders.

The government seems to be moving in this direction, but this long process should not prevent for the time being the continuation of capital accumulation, formal or informal, in residential real estate in the largest cities.

19 Some small towns experienced the bursting of a speculative bubble, but the authorities tried to stifle the information. This was notably the case of Wenzhou, a "small town" of 3 million inhabitants which experienced a serious real estate crisis in 2011. See Yufeng Chen, Ye Zhipeng, and Huang Guan (2018

  2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 Central gouvernement revenues Central government expanses Local governement expanses Local government revenues X100 million CNYis held by households as pure savings, developer inventories have multiplied in small and medium-sized cities affected by deindustrialization where the success story of the metropolises could not be replicated 5 .

Figure 2 .

 2 Figure 2. The mechanism of land value capture to finance urbanization

Figure 3 .Figure 4 .

 34 Figure 3. Residential vacancy in several cities (in percentage of the total housing stock)

  6 https://www.shangyexinzhi.com/article/2122354.html 7 This is especially true for Western funds, with Hong Kong and Singaporean operators being strongly favoured by their networks on the Chinese continent. See N. Aveline-Dubach (2020) « The Financialization of Real Estate in Megacities and its Variegated Trajectories in East Asia", In André Sorensen and Danielle Labbé (dir.), Handbook of Megacities and Megacity Regions, Edward Elgar, 394-409. https://halshs.archives-ouvertes.fr/halshs-02517518/document. 8 Numbeo database, https://www.numbeo.com/property-investment/in/Beijin 9 The hukou is a system for registering Chinese citizens according to their place of origin. Citizens of rural origin who live and work in the cities do not have access to the social benefits of urban hukou holders (education,

housing, health, pension...). However, the situation is changing with the introduction of a points-based hukou system. See Cinzia Losavio (2018), "Building wealth through a stratified inclusion: the point-based hukou system in Zhuhai". In: Aveline-Dubach, N. (ed.) Pathways of Sustainable Development Across China -The Cases of Hangzhou, Datong and Zhuhai, pp. 109-121. Venice Imago Editor, Venice. https://halshs.archivesouvertes.fr/halshs-02873431/document 10 https://journals.openedition.org/cybergeo/27655

  ), "The Financial Crisis in Wenzhou: An Unanticipated Consequence of China's 'Four Trillion Yuan Economic Stimulus Package'". China: An International Journal 16 (1): 152-173. 20 Natacha Aveline-Dubach (2020) « The Financialization of Real Estate in Megacities and its Variegated Trajectories in East Asia", In André Sorensen and Danielle Labbé (dir.), Handbook of Megacities and Megacity Regions, Edward Elgar, 394-409. https://halshs.archives-ouvertes.fr/halshs-02517518/document.

Theurillat, T., Lenzer, J. and H. Zhan (2018), 'The increasing financialization of China's urbanization', Issues and Studies ,

(4). http://www.fingeo.net/wordpress/wp-content/uploads/2017/07/WP6_Theurillat-et-al_FinancializationChina.pdf

Mingye Li (2017), Evolution of Chinese Ghost CitiesOpportunity for a Paradigm Shift? The Case of Changzhou, China Perspectives; 1,

69-78. https://journals.openedition.org/chinaperspectives/7209.

Adèle Esposito Andujar (2020) «Figures d'un urbanisme aspirationnel. Les villes secondaires d'Asie du Sud-Est à l'heure de la Belt and Road Iniative», Urbanités, https://www.revue-urbanites.fr/usea-esposito/

Buying a home is generally considered a prerequisite for marriage for young men. New homes are preferred to mark the transition to a new life.

https://asia.nikkei.com/Business/Business-trends/China-s-Vanke-pilots-build-to-lease-Beijing-housingproject.

14 https://www.letemps.ch/economie/dette-gouvernements-locaux-chine-une-bombe-

retardement 15 www.lincolninst.edu/publications/articles/chinas-property-tax-reform

These transactions take place in "urban villages", former agricultural enclaves caught up by urbanisation where land ownership is not under the control of the local government but of a collective of villagers. The latter circumvent the rules of land allocation by illegally selling flats to internal migrants at values 40-60% below the local market. See the entry "Xiaochanquan (China)" by Cinzia Losavio in the online encyclopedia "Global informality project", https://www.informality.com/wiki/index.php?title=Xiaochanquan_(China).

On the conditions for financing Chinese property developers, see Theurillat T. (2014), "The financing of urban production in China", working document, http://www2.unine.ch/socio/page-24681.html-24681.html