Guillaume Bellec

Mathieu Galtier

Romain Brette

Pierre Yger

Slow feature analysis

Keywords: unsupervised learning, plasticity, slow feature analysis

Introduction

A property of sensory cortices is the ability to extract and process invariants in the flow of information they are receiving. The higher the cognitive areas, the more complex the invariants [START_REF] R Quian Quiroga | Invariant visual representation by single neurons in the human brain[END_REF][START_REF] Wallis | Invariant face and object recognition in the visual system[END_REF]. Finding invariants, or equivalently seeking statistically relevant features in the inputs is an unsupervised learning task performed by the sensory cortices [START_REF] Marr | A theory for cerebral neocortex[END_REF]. In the context of object recognition, for example, given several images of a plane and several others of a car, learning to answer the question "is it a car?" when facing a new image is rather complex. To answer correctly, one must have identified an invariant belonging to the "car" category, common in all the images. By taking time into account, one could use the fact that two images appearing one after the other in a natural scene likely originate from the same object [START_REF] James J Dicarlo | How Does the Brain Solve Visual Object Recognition[END_REF]. This shows the importance for cortical areas to learn from temporal features, and how time can be used as an extra dimension to extract information.

To address more generally this question of time in unsupervised learning, one can start from the observation that sensory signals tend to vary at a fast time-scale compared to the source of these signals in the environment. For example, the physical location of a sound source varies more slowly than the acoustical pressure. This has led several authors to propose statistical learning algorithms that extract slowly varying features of sensory signals [START_REF] Wiskott | Slow Feature Analysis: unsupervised learning of invariances[END_REF][START_REF] Öldiák | Learning invariance from transformation sequences[END_REF][START_REF] Einhäuser | Learning viewpoint invariant object representations using a temporal coherence principle[END_REF][START_REF] Becker | Self-organizing neural network that discovers surfaces in random-dot stereograms[END_REF], which are also related to psychological theories according to which the basis of perception is the "invariant structure" of the sensory (or sensorimotor) signals [START_REF] Gibson | The Ecological Approach to Visual Perception[END_REF].

Slow Feature Analysis (SFA) [START_REF] Wiskott | Slow Feature Analysis: unsupervised learning of invariances[END_REF][START_REF] Wiskott | Is slowness a learning principle of the visual cortex?[END_REF] is a simple and yet very successful unsupervised algorithm to capture invariant features from input streams. Literature on SFA provides arguments that the brain might process information in a similar way as SFA. In the primary visual cortex V1, the algorithm has been used to explain the emergence of simple and complex cells, and the observed diversity of their properties [START_REF] Berkes | Slow Feature Analysis yields a rich repertoire of complex cell properties[END_REF][START_REF] Kayser | Extracting Slow Subspaces from Natural Videos Leads to Complex Cells[END_REF][START_REF] Dähne | Slow Feature Analysis on retinal waves leads to V1 complex cells[END_REF]. In the somatosensory system, it has also been used to show the emergence of place and head-direction cells [START_REF] Franzius | Slowness and sparseness lead to place, head-direction, and spatialview cells[END_REF]. SFA extracts continuous causes of the inputs [START_REF] Turner | A maximum-likelihood interpretation for Slow Feature Analysis[END_REF] and such stochastic hypotheses also seem to be compatible with known phenomenon of neural adaptation [START_REF] Stevenson | Sensory adaptation and short term plasticity as bayesian correction for a changing brain[END_REF]. In this paper, we will show how those principles can be applied to inputs where time is intrinsically entangled to the signals themselves, i.e in the auditory pathway.

If the brain is able to achieve such representations through unsupervised learning, how is it implemented with spiking neurons? While it has already been shown [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF] that Slow Feature Analysis (SFA) can be linked to the trace learning rule [START_REF] Öldiák | Learning invariance from transformation sequences[END_REF], and that an implementation of such a trace rule can be theoretically compatible with known experimental facts on synaptic plasticity such as Spike-Timing Dependent Plasticity (STDP) [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF][START_REF] Markram | Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs[END_REF][START_REF] Gerstner | A neuronal learning rule for sub-millisecond temporal coding[END_REF], a simulation of rate-based or spiking neurons implementing SFA is still missing. In this work, we show mathematically in a similar way as in [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF] and with simulations of rate-based and spiking neurons how an optimization criterion derived from slowness can be implemented by neurons, and we apply this in the context of auditory signals.

While the derivation is similar to [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF], the conclusions are slightly different. As in their results, we show that synapses can be governed by a biologically plausible online learning rule capturing the slowest time-varying component of the inputs. Yet substantial extensions of their work are added since we found through theory and comparative simulations that classical asymmetric STDP [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF] and trace learning [START_REF] Öldiák | Learning invariance from transformation sequences[END_REF] are not sufficient to implement SFA, when sharp postsynaptic potentials are considered, while a plasticity kernel shaped as a smoothed second derivative is.

Material and methods

Notations

In all the following x(t) ∈ R n denotes the input signal and s is the output such that ∀t, s(t) := w T ν(t), w = [w 1 , .., w n ] being the incoming synaptic weights. To simplify the theory all signals are considered null outside a finite support [0, T]. The convolution between signals x and y is written x * y, the integration over time is written x, and cross-correlation is written x, y := T 0 x(t)y(t)dt. We write g τ (t) for the temporal filter given by g τ (t) := 1 τ H(t)e -t τ , where H is the Heaviside function and τ is a time constant. From the mathematical framework detailed in [START_REF] Galtier | Multiscale analysis of slow-fast neuronal learning models with noise[END_REF] we interpret convolution between temporal signals as simple matrix multiplication. The transposition notation is extended to the continuous representation of filters by noting g T τ (t) = g τ (-t) because the property: x * g, y = x, y * g T will be widely used in our paper. The derivative operator at order (r) is also written as a convolution with d (r) such that x * d (r) := d r x dt r . δ t 0 is a Dirac distribution centred at time t 0 .

Toy example

To test the convergence of SFA, we use the toy example defined in [START_REF] Wiskott | Slow Feature Analysis: unsupervised learning of invariances[END_REF]. Five inputs x i (t) are built out of non linear combinations of two sinusoids:

               x 1 (t) = sin(2π f 0 t) + α cos(2π11 f 0 t) 2 x 2 (t) = cos(2π11 f 0 t) x 3 (t) = x 1 (t) 2 x 4 (t) = x 1 (t)x 2 (t) x 5 (t) = x 2 (t) 2 (1)
In all the manuscript, we consider f 0 = 1Hz.

The spiking neuron model

r i := θ i * φ
Note that r i is in Hz. For simplicity, all input rates ν i are scaled to have a time average of r and a standard deviation r s so that ν i := r + r s x i . We assume as in [START_REF] Kempter | Hebbian learning and spiking neurons[END_REF][START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF][START_REF] Eugene | Relating STDP to BCM[END_REF] for modelling STDP that the post-synaptic neuron generates a Poisson spike train θ out of rate ν out (t) that is a linear function of the post synaptic potential (PSP) produced by given spike trains θ i . In other words, we have a frequency offset ν 0 , a constant κ, synaptic weights w i and a normalized PSP filter ξ (ξ = 1) that defines the neuron dynamic as

ν out (t) = ν 0 + κ ∑ i w i θ i * ξ (2)
The plausibility of this assumption is discussed in [START_REF] Kempter | Hebbian learning and spiking neurons[END_REF][START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF], and mathematical details are given in the Appendix. For numerical simulations, we chose r = 100 Hz, r s = 80 Hz, ν 0 = 100 Hz and κ = 0.0625. In Figure 6 τ rate = 10 ms and 2 ms in Figure 9. The PSP filter used in Figures 6 and 9 is a decreasing exponential of time constant 1 ms. Simulations have been performed with the Brian simulator [START_REF] Goodman | Brian: A Simulator for Spiking Neural Networks in Python[END_REF], using a time step δt = 0.1 ms.

Time scales assumptions

We assume that if δt is the discrete time step of the simulation, τ rate the width of the spike counting filter, τ signal the typical time scale of variation of the useful information hidden in the input, τ w the time constant of synaptic updates, we have

In the following, τ STDP the width of the plasticity kernel is equal to τ rate . Inequality (3.b) is needed so that an online implementation converges towards the same solution as a batch implementation and (3.a) guarantees that information encoded into the rate does not vary faster than the time required to compute the rate.

Calculation of Ω SFA The exact shape of Ω SFA is calculated starting from the expression [START_REF] Becker | Self-organizing neural network that discovers surfaces in random-dot stereograms[END_REF] . And using φ := g τ rate * g τ rate , we have

Ω SFA := φ * d (2) * φ T When reversing time d (1) becomes -d T (1) , thus d (2) = -d (1) * d T (
Ω SFA = -g τ rate * g τ rate * d (1) * d T (1) * g T τ rate * g T τ rate
The convolution is commutative so we can group the terms using the smoothed derivative filter

d (1,τ) := g τ * d (1) * g T τ .
This leads to

Ω SFA = -d (1,τ rate ) * d T (1,τ rate )
Proof is given in [START_REF] Mathieu | A biological gradient descent for prediction through a combination of stdp and homeostatic plasticity[END_REF] that

d (1,τ) = g T τ -g τ 2 .
So when developing the expression of d (1,τ rate ) we have , Λ Hebbian+ = δ 0 and Λ Hebbian-= -δ 0 . Those three filters are implementing a classical STDP like plasticity kernel (asymetric), a purely Hebbian plasticity kernel, and an anti-Hebbian one, respectively. They lead to the three control plasticity kernels

Ω SFA = - 1 4 (g
Ω Classic := d (1,τ) = H(t)e (-t τ ) -H(-t)e ( t τ )
2τ

Ω Hebbian+ := σ τ = e -|t| τ 2τ Ω Hebbian-:= -σ τ = - e -|t| τ 2τ
Measure of slowness For each set of parameters and STDP kernel, the batch algorithm is run 20 times on the toy example. Slowness of trial i is measured by the correlation coefficient CC i between output s i and the optimal solution sin(2π f 0 t), then the CC i are squared and averaged geometrically such that we can quantify the similarity to the optimal solution in Figure 8 as

CC = Π 20 k=1 CC 2/20
k . Note that if τ STDP = 0 ms, rate is ill defined and we replaced the spike train θ i by its ideal expected value.

Application to audio processing In Figure 9-11, we used recording sounds x 0 sampled at f s . Inputs to SFA algorithm are built by constructing delay lines of N channels x i (t) delaying x 0 by delay iτ delay such that

x i (t) = x 0 (t -iτ delay ) (5)
To be more precise, in Figure 9 f s = 11 kHz, N = 256, τ delay = 0.81 ms. To study the influence of parameters, the same simulation is run in Figure 10 with τ delay = 8.8 ms. In Figure 11 to capture both pitch and tempo we use N = 512 delays with τ delay = 1.5 ms.

Results

A number of models have been proposed to explain how unsupervised learning of invariants could be achieved by spiking neurons, either with a normative approach [START_REF] Toyoizumi | Optimality Model of Unsupervised Spike-Timing-Dependent Plasticity: Synaptic Memory and Weight Distribution[END_REF][START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF][START_REF] Rao | Spike-timingdependent Hebbian plasticity as temporal difference learning[END_REF] or with a phenomenological one [START_REF] Clopath | Connectivity reflects coding: a model of voltage-based STDP with homeostasis[END_REF][START_REF] Yger | The Convallis rule for unsupervised learning in cortical networks[END_REF][START_REF] Yger | Stable learning in stochastic network states[END_REF]. Slow Feature Analysis (SFA) [START_REF] Wiskott | Slow Feature Analysis: unsupervised learning of invariances[END_REF][START_REF] Wiskott | Is slowness a learning principle of the visual cortex?[END_REF] is a simple and yet very successful unsupervised algorithm to capture invariant features. In Figure 1A

x(t)

Time Fig. 1 Illustration of the SFA algorithm and Spike-Timing-Dependent plasticity. A If x are the sensory inputs, SFA learns a linear transformation W that maximizes the slowness of the output s. s(t) depends only on x(t) and not on its past. B Illustration of the Spike-Timing-Dependent Plasticity protocol, as observed in vitro in [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF][START_REF] Markram | Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs[END_REF]. The synaptic weight is modified as a function of the temporal difference ∆t between pre and post synaptic spikes.

of the SFA algorithm is illustrated. If x i (t) are the inputs, SFA will decompose it onto a basis of signals s i (t) sorted by their slowness. In this paper, we focus on how neurons can implement such an algorithm. In line with seminal ideas of Hebbian learning [START_REF] Do O Hebb | The organization of behavior: a neuropsychological theory[END_REF], and since its first discovery [START_REF] Bliss | Long-lasting potentiation of synaptic transmission in the dentate area of the anaesthetized rabbit following stimulation of the preforant path[END_REF], Spike-Timing dependent plasticity (STDP) appears to be a good candidate to support such a dynamical extraction of invariant by neurons, as a widely observed phenomenon taking place in a large number of systems (see [START_REF] Dan | Spike Timing-Dependent Plasticity of Neural Circuits[END_REF] for review). As illustrated in Figure 1B, it relies on the fact that fine pre-post interactions [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF][START_REF] Markram | Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs[END_REF][START_REF] Gerstner | A neuronal learning rule for sub-millisecond temporal coding[END_REF] can modify the strength of a single synaptic weight. It has already been stated that some STDP kernel (Figure 1B) can be used to implement SFA [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF]. But the question of the most efficient plasticity kernel implementing it remains. Here, we compute analytically what should be the shape of a plasticity kernel implementing SFA and provide simulations of its implementation to compare with kernels suggested in [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF].

Implementing SFA with gradient methods

We consider the theoretical formulation of SFA as an optimisation problem. Originally [START_REF] Wiskott | Slow Feature Analysis: unsupervised learning of invariances[END_REF] the outputs of SFA are defined as multi-dimensional and decorrelated with each other. However as find-ing slow components and decorrelating them are two distinct tasks, here we focus only on the first (for the decorrelation task see [START_REF] Pozzorini | Temporal whitening by power-law adaptation in neocortical neurons[END_REF][START_REF] Dan | Efficient coding of natural scenes in the lateral geniculate nucleus: experimental test of a computational theory[END_REF]). Therefore, we assume as in [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF] that the inputs x are normalised and decorrelated so that x, x = Id N . Note that in all simulations inputs are sphered during a pre-processing stage to satisfy this constraint and this is not implemented biologically in this paper (see Discussion). Given those N input signals x i (t) with zero mean and unit variance, the aim of SFA is to find a set of weights w maximising the slowness of the output s(t) = w T x(t). This problem is formalised in [START_REF] Wiskott | Slow Feature Analysis: unsupervised learning of invariances[END_REF] as min w ds dt 2 subject to:

s 2 = 1 ∀t, s(t) = w T x(t) (6)
To illustrate the algorithm, we define a toy example (see Methods), and as can be seen in Figure 2, SFA isolates in the output s the slow component of the inputs. In all the following, this simple example will be used to compare with SFA.

SFA can be re-written so that we can derive a projected gradient algorithm to solve the problem. Because the covariance of the inputs is identity, normalising the weights is equivalent to normalising the output variance. After integrating by part the squared derivative of s, the cost function becomes -d 2 s dt 2 , s1 . Because s = w T x, and using the notation of the second order derivative operator as a convolution with d (2) (see Methods), we

have d 2 s dt 2 = w T x * d (2)
. Assuming weights are constant between iterations of a batch algorithm, SFA can be reduced to finding w such that max

w w T x * d (2) , x w subject to: ||w|| = 1 (7)
As far as the objective is quadratic in w, the gradient is easily written as 2 x * d (2) , x w. Beside normalization at each time step, the iteration of the batch algorithm is given by

w k -w k-1 = ∆w k ∝ x * d (2) , s (8)
A key issue of the batch algorithm is that the time flow is not realistic (see Figure 3A). In the batch algorithm, at each iteration the exact same stimulus is presented and the weights are updated based on statistics of the whole sequence. Instead,

(t) ∝ d 2 x i dt 2 (t)s(t).
One can note that the only difference as compared to Oja's implementation of PCA (Principal Component Analysis) ẇi (t) ∝ x i (t)s(t) [START_REF] Oja | Simplified neuron model as a principal component analyzer[END_REF] is the use of the second derivative. The same result could have been obtained by keeping one derivative on each side ( ẇi ∝ ẋi ṡ) but grouping the two make sense when dealing with STDP in the following sections.

In Figure 4, we compare the different implementations of the batch and the online algorithms on the toy example defined in Figure 2. As compared to the ground truth given by the standard SFA algorithm (Figure 4A), the batch (Figure 4B) and the online (Figure 4C) implementations all converge toward the same solution. Moreover, in the case of the online implementation, the dynamics of the weights (Figure 4D) shows that the convergence is achieved in a few seconds of simulated time.

To show that recovering the slow sinusoid in this example is due to the use of the second derivative d [START_REF] Berkes | Slow Feature Analysis yields a rich repertoire of complex cell properties[END_REF] in Equation ( 7), we compared different learning rules on the exact same task (see A Slowest component given by SFA as in [START_REF] Wiskott | Slow Feature Analysis: unsupervised learning of invariances[END_REF] for the toy example described in Figure 2. B Same but for the batch implementation. (Phases of the oscillations are opposed, but the two solutions are equally optimal). C Same for the online form: it converges to the same solution in a few seconds. D Evolution of the five corresponding weights as a function of time with the online implementation. Convergence is reached in a few seconds.

Methods) by replacing d (2) in the batch and online algorithm by Λ Classic := d (1) , Λ Hebbian+ := δ 0 and Λ Hebbian-:= -δ 0 . Classic stands for the fact that this is at the origin of the canonical kernel for plasticity, as described in Figure 1B [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF][START_REF] Markram | Regulation of synaptic efficacy by coincidence of postsynaptic APs and EPSPs[END_REF][START_REF] Gerstner | A neuronal learning rule for sub-millisecond temporal coding[END_REF]. Note that the Hebbian learning rule using Λ Hebbian+ , if normalized properly, is known to recover Oja's neural implementation of PCA [START_REF] Oja | Simplified neuron model as a principal component analyzer[END_REF], defined by ẇi ∝ x i * δ 0 (t)s(t) = x i (t)s(t). The left column of Figure 5 shows the pattern of the output after convergence of the batch algorithm, while the right column shows convergence of the online implementation. Only Λ SFA = d (2) is able to recover the slow sinusoid.

Implementation by spiking neurons

To propose an implementation of SFA with biological neurons, we need to take into account the spiking property of real neurons. Input x i is scaled between 20 and 180 Hz to be encoded by pre-synaptic firing rates ν i . It generates the Poisson spike train θ i and we can define a rate estimator r i := θ i * φ rate , built by counting spikes over a short time bin (see Methods). We assumed, such as in [START_REF] Kempter | Hebbian learning and spiking neurons[END_REF][START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF][START_REF] Eugene | Relating STDP to BCM[END_REF], that the post-synaptic neuron generates a Poisson spike train θ out of rate ν out (t) linearly related to the post synaptic potential (see Methods). The rate estimator r out (t) is defined as for the pre-synaptic spike trains.

In the following we show that for a kernel shaped as a smoothed second derivative, the ex- pected weight modification due to STDP implements the gradient ascent of Equation [START_REF] Dan | Spike Timing-Dependent Plasticity of Neural Circuits[END_REF]. Let t in i,n (resp. t out m ) be the spike times of input neuron i (resp. output neuron). If the filter that defines the STDP plasticity kernel is called Ω, and assuming synaptic changes are occurring on a slow timescale τ w compared to the typical interspike interval, we have, if is a small learning rate

t t-τ w ẇi (t)dt = ∑ t in ,t out ∈[t-τ w ,t] Ω t out -t in i = t t-τ w θ i * Ω(t )θ out (t )dt Such that, on average ẇi = θ i * Ω, θ out (9)
We assume the existence of a kernel Λ so that Ω can be seen as a smoothed version of it, i.e. Ω = φ rate * Λ * φ T rate . In particular, the width of such a plasticity kernel Ω is controlled by φ rate . If Λ = δ 0 , then the width of the plasticity kernel is 2τ rate . In all the following, to be consistent, we define τ STDP as controlling the width of this STDP kernel (τ STDP = τ rate ). From Equation ( 9), we have [START_REF] Einhäuser | Learning viewpoint invariant object representations using a temporal coherence principle[END_REF] We see that for Λ = d (2) , the average effect of our learning rule built with an STDP-like plasticity kernel has a similar effect as the Hebbian learning rule ẇi ∝ d 2 x i dt 2 s that implements SFA. To generalize the initial batch and online algorithms without spikes 2 we use

ẇi = θ i * φ rate * Λ * φ T rate , θ out = θ i * φ rate * Λ, θ out * φ rate = r i * Λ, r out
Ω SFA = φ rate * d (2) * φ T rate
It implies that the plasticity kernel Ω SFA has to be shaped as a smoothed second derivative, and the exact shape of Ω SFA depends only on φ rate . The calculation that exhibits the literal expression of the kernel Ω SFA is detailed in the methods, and it leads to the plasticity kernel below

Ω SFA (t) = 1 2τ STDP e -|t| τ STDP |t| τ STDP -1 (11)
This demonstration gives a good intuition of the reason why this plasticity kernel Ω SFA (t) implements SFA. However, to see the full mathematical derivation, the reader should refer to the Appendix, where the stochastic aspect of the input and output spiking activities are taken into account, as in [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF] using results from [START_REF] Kempter | Hebbian learning and spiking neurons[END_REF]. The outcome is that to implement SFA, the expected weight modification due to STDP has to be proportional to x i * Ω SFA , s . This is verified by computing the stochastic average of Equation [START_REF] Einhäuser | Learning viewpoint invariant object representations using a temporal coherence principle[END_REF]. Beside terms that can be neglected (see Appendix or [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF]) the expected weight change is in fact proportional to x i * Ω SFA * ξ T , s , where ξ is the post synaptic potential filter. But if the time constant 2 The non-spiking case is recovered if τ rate → 0 and θ i (t) is replaced by the expected value being the rate ν of the inputs. of ξ is small compared to τ rate , the filtering performed by ξ T can be neglected and STDP using Ω SFA implements SFA. In simulations we used time constants of the order of 1 ms for both ξ and Ω to match biological data and the results are not impaired.

In Figure 6, we test the implementation of SFA with spiking neurons. Figure 6A illustrates the spiking property of the neuron. Variations of the output rate estimator r out over time are represented in Figure 6B. After convergence, the shape of the sinusoid is approximately being reproduced. The noise is large because the input neurons are spiking at reasonably low rates (between 20-180 Hz, see Methods). To confirm that the slow sinusoid is well captured, we use the periodicity of the stimulus to average multiple chunks of two seconds length. When averaged over all periods after convergence, the averaged pattern can be represented as in Figure 6C). [START_REF] Oja | Simplified neuron model as a principal component analyzer[END_REF] when τ STDP → 0 ms and equivalent to the trace learning [START_REF] Öldiák | Learning invariance from transformation sequences[END_REF] in general. D Anti-Hebbian kernel.

Figure 7. In particular Ω Hebbian+ is equivalent to trace learning [START_REF] Öldiák | Learning invariance from transformation sequences[END_REF] according to our theory and Ω Classic , measured with biological neurons in vitro [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF] is mostly used in the literature. We show in Figure 8 that only Ω SFA implements SFA for any choice of parameters.

The parameter α controlling the amplitude of the quickly varying components in the toy example is varied for different values of the STDP time constant τ STDP . The influence of α is shown in Figure 8A, while Figure 8B shows the modulation of the plasticity kernel as a function of τ STDP = τ rate . Figure 8C shows which STDP kernel is capable of recovering the slow sinusoid with an extreme case of a width of τ STDP = 0 ms (left), τ STDP = 10 ms (middle) and α = 1, τ STDP being varied (right). Our implementation of SFA recovers the slow solution in each case (see Methods). For large realistic STDP kernel width, the Hebbian kernel only recovers the slow sinusoid for α < 1000 and the anti-Hebbian for α > 1000. Instead Ω SFA finds the expected solution for any α. When varying τ STDP , there is only one case where Λ SFA does not recover the slow sinusoid: when the plasticity kernel is larger than the fastest oscillations (τ STDP = 100 ms, right panel, last row). But because the assumption of the time scales (3.b, see Methods) is not respected any more, we conclude that the slowest component is robustly extracted by a plasticity kernel shaped by Ω SFA only.

Application to audio recordings

Having shown the validity and the robustness of our implementations on a simple example, we can now use more complex inputs. Because in audio time is crucial and because SFA has mostly been applied to visual stimuli [START_REF] Berkes | Slow Feature Analysis yields a rich repertoire of complex cell properties[END_REF][START_REF] Kayser | Extracting Slow Subspaces from Natural Videos Leads to Complex Cells[END_REF][START_REF] Dähne | Slow Feature Analysis on retinal waves leads to V1 complex cells[END_REF], music and speech recordings are good candidates to try our implementations of SFA.

To begin with, we apply SFA on a guitar signal x 0 (15s of a country music), sampled at f s = 11 kHz with our model of a plastic spiking neuron (see Figure 9A for a schematic representation). This one dimensional audio recording is expanded to dimension 64 by delay lines (see Methods), providing multiple delayed versions of x 0 , each of them delayed by 0.81 ms compared to the previous one (see Figure 9A). From each channel x i an inhomogeneous time-varying Poisson spike train θ i is generated with rates varying between 20 and 180 Hz (see Methods). Using the STDP kernel Ω SFA we see on Figure 9B the convergence of the weights from all those channels when looping over that country guitar recording. Figure 9C displays the final weights assigned to each delay, i.e. the learnt filter given by SFA. We can see that such an optimal filter is a sinusoid at a particular frequency, which corresponds to the dominant frequency of the recording (97Hz being the frequency of G2 the note mostly played by the guitarist).

By using delay lines with larger delays (8.8 ms) between each inputs (see Figure 10A), we can see on Figure 10B the results obtained with the standard SFA algorithm and our implementation. We can see in frequency domain (Figure 10C) that the solutions peak at 3 Hz. Filters, again, are sinusoids (Figure 10D) tuned to about 3 Hz (Figure 10E), i.e. to 180 beats per minute, the tempo of the music piece. Interestingly, the visualisation of the filters in time domain in Figure 10D seems to show that the batch implementation produced a more regular solution than the classic SFA algorithm.

Assuming the length of the delays in the delay lines can impact the filters obtained by SFA algorithm, we ask what is the influence of the parameter τ STDP , the width of the plasticity kernel. In Figure 11, we used the same guitar recording as previously expanded with 512 delays of 1.5 ms, to get a broad range of possible filters. The plasticity kernel has a variable width τ STDP varying between 0.1 ms to about 30 ms. Figure 11A,B

Learnt lter

Fig. 10 Extracting the tempo of a guitar signal. A Schematic construction of the inputs, as in Figure 9. A single guitar recording x 0 of country music is delayed 256 times into channels x i (t) with delays of multiple of τ delay = 8.8 ms (delay line), and fed into the SFA algorithm. B The audio recording x 0 (blue), compared with the two outputs of SFA (filtered versions of input). In black, the output is obtained with our gradient based implementation of SFA, and in red, this is the original algorithm (ground-truth). C Comparison of the power spectrum for the three signals, showing an increase in energy at the tempo frequency (3 Hz). D Weights after learning for all input channels with the two implementations of SFA. E The power spectrum of the two learnt filters. They are tuned to 3 Hz, the tempo of the recording.

is, the lower the captured frequencies are. Figure 11D represents, in a more exhaustive manner, the frequency responses as a function of τ STDP . Only few frequencies have high energy in the learnt frequency response, and the larger the kernel is, the lower the tuning frequencies of the filter are.

For sharp STDP kernel (τ rate ≈ δt) the pitch G2 is captured. With a slightly larger one the filter is tuned at the pitch E2 (another bass note played by the musician). For wide STDP kernel, the tempo is captured even with small delays between the inputs. In between those note pitches and the tempo, the captured frequencies are more complex to interpret.

Discussion

In this work, we developed a mathematical framework showing how SFA could be implemented by models of neurons. We first designed an online learning rule to extract the slowest components of their inputs, and then adapted it in the context of spiking neurons. It pro-vides a neural implementation of SFA compatible with the experimental observation of Spike-Time-Dependent-Plasticity [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF][START_REF] Gerstner | A neuronal learning rule for sub-millisecond temporal coding[END_REF] and rate homoeostasis. By implementing a gradient descent on the cost function of SFA, we found the precise plasticity kernel that fits a particular function. We showed through simulations that the STDP kernel has to be shaped as a reversed Mexican hat (second derivative) to learn the optimally slow features. When applied to auditory signals, simulations showed that our gradient implementation of SFA allows a single neuron to extract the tempo of its inputs.

The compatibility of SFA and STDP was already demonstrated in [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF]. With continuous and deterministic inputs they characterized that with any filter Λ having a Fourier transform of the formf 2 ( f is the frequency) truncated to be zero at frequencies higher than f max , a gradient descent of the form ∆w k ∝ x * Λ, s implements SFA when the input information does not contain frequencies higher than f max . Because the second derivative is a convolution by a filter beingf 2 in Fourier domain, our analysis is strictly equivalent but does not require to limit the input frequencies below f max .

In this work, various kernels have been tested as controls, and among those, the classical STDP kernel Ω Classic and the trace learning equivalent to a plasticity kernel using Ω Hebbian+ . In agreement and as suggested in [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF], we found that the trace learning rule is a good approximation of SFA (see performances on Figure 8C) for simulations on several toy examples. However, this control kernel does not extract the slow feature any more when a lot of energy is contained in the high frequencies. Control simulations with Ω classic and the insight of our implementation with a second order derivative allow to raise the problem that the classical asymmetric STDP described in [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF] is not suited for SFA. The classical asymmetric STDP kernel was supported in [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF], and the divergence comes from the fact that they used PSP time constants much larger than τ rate when we did the opposite. In the limit of large ξ, we also find 3 an asymmetric kernel but it is a smoothed third or-der derivative which is not as trivial as the classic STDP.

In all the paper, the expression of the SFA problem is simplified assuming inputs have been already pre-processed and decorrelated. Various theories may explain how inputs onto sensory cortices can be decorrelated [START_REF] Pozzorini | Temporal whitening by power-law adaptation in neocortical neurons[END_REF][START_REF] Dan | Efficient coding of natural scenes in the lateral geniculate nucleus: experimental test of a computational theory[END_REF], but a plausible implementation of such a process with spiking neurons is beyond the scope of this manuscript. In addition, a clear restriction here is that in the generic formulation of SFA, a decorrelation process at the output stage is required to provide a rich population of non-redundant SFA neurons. However, we believe that such a decomposition could be mediated with lateral inhibition and competition between several output neurons [START_REF] Michael S Falconbridge | A simple Hebbian/anti-Hebbian network learns the sparse, independent components of natural images[END_REF].

Appendix

Filtered auto-correlation for a Poisson spike train

The main difference between the implementation of SFA with a spiking neuron and with a linear-rate-based neuron raises from the fact that the expected value of the auto-correlation of a Poisson spike train is not the auto-correlation of the rates; an additional term appears. This is justified in [START_REF] Kempter | Hebbian learning and spiking neurons[END_REF] and used in [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF]. Formally, for any Poisson spike train θ of rate ν, time t, given two filters f and g, we have if E(x) stands for the stochastic average over multiple realizations of a random variable x

E [ θ * f , θ * g ] = ν * f , ν * g + ν f g (12)
Equation ( 12) is obtained because we have ∀u, vE [θ(u)θ(v)] = ν(u)ν(v) + δ u (v)ν(u) because spiking at time u and v are dependent only if u equals v for which the variance of θ(u) is equal to its mean ν(u). Switching filtering integration and ensemble integration we have E [θ * f (t)θ * g(t)] = u v g(tu) f (tv)E [θ(u)θ(v)] dvdu. Inserting the previous equality leads to

E [θ * f (t)θ(t)] = ν * f (t)ν * g(t) + ν * ( f g)(t)
and integration over time using the fact that ν * ( f g) = ν f g gives Equation [START_REF] Öldiák | Learning invariance from transformation sequences[END_REF].

Average effect of STDP for spiking neurons From Equation (9) we can write the net effect of STDP ẇ = θ i * Ω, θ out To compute the stochastic average of STDP over multiple realizations of the spike trains, we use

E [ ẇ ] = E [E [ ẇ |θ i ]].
In other words, for a fixed set of input spike trains θ i we first integrate over all the possible output spike trains, we insert the linear dynamics of the neuron in Equation ( 2) and only terms constant in θ remain

E [ ẇi |θ i ] = E [ θ * Ω, θ out |θ i ]
= κ ∑ j w j θ i * Ω, θ j * ξ + ν 0 θ i * Ω Averaging now over all the possible input spike trains, θ i and θ j being independent for i = j, we can replace θ i and θ j by their expected value ν i and ν j . For the remaining term i = j we use Equation [START_REF] Öldiák | Learning invariance from transformation sequences[END_REF], such that we have

E [ ẇi ] = κ n ∑ j=1 w j ν i * Ω, ν j * ξ + κw i rΩξ + ν 0 rΩ (13)
The third term can be neglected as in [START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF], in our case Ω is null anyway. The second remaining term is proportional to w, and this is equivalent to simply adding in the cost function a term proportional to w 2 which would be made obsolete by the weight normalization. What remains for the effect of STDP on average is therefore proportional to ∑ j w j ν i * Ω, ν j * ξ which can be rewritten as ν i * Ω * ξ T , ∑ j w j ν j . Finally, to rigorously show that this is equivalent to the gradient implementation of SFA with continuous neurons obtained in Equation 8, we replace ν by r + r s x in Equation 13and we obtain E [ ẇi ] = κr 2 s x i * Ω * ξ T , ∑ j w j x j + κr 2 s x i * Ω * ξ T , rw sum + κr 2 s rΩ ξ T , ∑ j w j ν j

The second term is null because it is equivalent to integrating a filtered version of a derivative of x over time and x is null at infinity, and the third term is null because Ω is null. What remains is therefore only

E [ ẇi ] = κr 2 s x i * Ω * ξ T , s (14)

2 Fig. 2 Fig. 3

 223 Fig.2Illustration of SFA. Simple toy example for the SFA algorithm, as described in[START_REF] Wiskott | Slow Feature Analysis: unsupervised learning of invariances[END_REF]. Five inputs (in blue) are nonlinear combinations of a slow and a fast sinusoids (see Methods). Output of the SFA algorithm is shown in red. Whitening of the inputs is not shown for clarity, but we have x, x = Id N

Fig. 4

 4 Fig.[START_REF] Bliss | Long-lasting potentiation of synaptic transmission in the dentate area of the anaesthetized rabbit following stimulation of the preforant path[END_REF] Convergence of the algorithms on the toy example. A Slowest component given by SFA as in[START_REF] Wiskott | Slow Feature Analysis: unsupervised learning of invariances[END_REF] for the toy example described in Figure2. B Same but for the batch implementation. (Phases of the oscillations are opposed, but the two solutions are equally optimal). C Same for the online form: it converges to the same solution in a few seconds. D Evolution of the five corresponding weights as a function of time with the online implementation. Convergence is reached in a few seconds.

Fig. 5

 5 Fig.[START_REF] Clopath | Connectivity reflects coding: a model of voltage-based STDP with homeostasis[END_REF] Comparison with other control learning rules. Four different rules are compared (rows) on the toy example described in Figure2. Left column shows the output pattern after convergence of the batch algorithm, right shows the results obtained with the online algorithm. A SFA kernel as a second derivative. B Plasticity kernel based on first order derivative. C Hebbian kernel rule implementing PCA[START_REF] Oja | Simplified neuron model as a principal component analyzer[END_REF]. D Anti-Hebbian kernel.

Fig. 6

 6 Fig.[START_REF] Dähne | Slow Feature Analysis on retinal waves leads to V1 complex cells[END_REF] Generalisation of the online implementation to ratebased neurons. A We consider r out as the output of the neuron to represent the slow output. B Implementation of SFA for the toy example described in Figure2with five inputs being non-homogeneous Poisson processes of rates x 1 , . . . x 5 re-scaled to vary from 20 to 180 Hz. Evolution of r out during learning. C Average over multiple chunks of 2 seconds showing the mean pattern extracted by the algorithm. D Evolution of the five synaptic weights during learning, to show that convergence is achieved.

Fig. 7

 7 Fig.[START_REF] Dan | Efficient coding of natural scenes in the lateral geniculate nucleus: experimental test of a computational theory[END_REF] Representation of the STDP kernels implementing the control learning rules for plasticity (see Methods for literal expressions). A SFA kernel as a second derivative. B Plasticity kernel based on first order derivative. C Hebbian kernel rule implementing PCA[START_REF] Oja | Simplified neuron model as a principal component analyzer[END_REF] when τ STDP → 0 ms and equivalent to the trace learning[START_REF] Öldiák | Learning invariance from transformation sequences[END_REF] in general. D Anti-Hebbian kernel.

BFig. 8 Fig. 9

 89 Fig. 8 Robustness to noise and to parameters of the plasticity kernels. A Influence of α, balancing the amount of high frequency information in the input (See Methods on the toy example for more details). B Influence of τ STDP , the time constant defining the width of the STDP kernel. C After fixing one parameter (Left τ STDP = 0 ms, Middle τ STDP = 10 ms, and Right α = 1) the other parameter (either τ STDP or α is varying on the horizontal axis). In each case the four STDP kernels introduced in Figure 7 are compared on the vertical axis. For each point of the graphs 20 batch simulations are launched with different random weight initialisation, and slowness is measured as a correlation with the slow sinusoid averaged over trials (see Methods).

Fig. 11

 11 Fig. 11 Influence of the width of the plasticity kernel τ STDP . Simulations are run on the guitar example (see Figure 10) with 512 delays and τ delay = 1.5 ms. Over the trials the dependency to the parameter τ STDP is explored with the batch implementation. A, B, C Three filters learnt with τ STDP at different orders of magnitude, from top to bottom τ STDP = 10 ms, τ STDP = 1 ms, τ STDP = 0.1 ms respectively. The larger τ STDP the lower the captured periodicity. D Power spectrum of the learnt filters, as a function of the width of the plasticity kernel, τ STDP , varied on the vertical axis. Each row is the frequency response of the learnt filter by SFA, color-coded. Dashed lines annotate the different identified peaks of the frequency response. In green, E2 and G2 are the two bass notes mostly played by the musician in the recording.

Construction of the plasticity kernels In

  order to test the performance of the plasticity kernel Ω SFA we compare it to control learning rules. To do so, knowing that Ω SFA is defined as φ * Λ SFA * φ T , with Λ SFA = d (2) , we change that filter Λ SFA by Λ Classic = d[START_REF] Becker | Self-organizing neural network that discovers surfaces in random-dot stereograms[END_REF]

		T τ rate -g τ rate ) * (g τ rate -g T τ rate )
	Developing again, only terms in g * g, g * g T or their transpose remain. When inserting the lit-
	eral expressions g τ * g T τ = 1 2τ e -|t| τ and g τ * g τ = H(t) t τ 2 e -t τ , we have
	Ω SFA (t) =	1 τ rate	e	-|t| τ rate	|t| τ rate	-1	(4)

  , the generic principle

	A		s(t)
	x 1
	x 2		S=WX	s (t) 1
	. .			s (t) 2
	x N .			3 s (t)
	Time
	B
	Spike Timing Dependent Plasticity
	t	t +	∆t
	Pre		Post
			∆t

The integration by part gives: ṡ, ṡ = [ ṡs] ∞ -∞s, s and the first term is null because the inputs are assumed to be null at infinity.

As done in[START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF] we searched for the STDP kernel that implements SFA by verifying d (2) = Ω * ξ T with ξ very large. In this limit, to invert the convolution by ξ we can perform a derivation (see[START_REF] Sprekeler | Slowness: An objective for Spike-Timing-Dependent Plasticity[END_REF]) so that the solution for Ω is a third order derivative.


Amine Aribi
email: amine.aribi@lmpt.univ-tours.fr

Ahmad El Soufi
email: ahmad.elsoufi@lmpt.univ-tours.fr

THE FIRST POSITIVE EIGENVALUE OF THE SUB-LAPLACIAN ON CR SPHERES

Keywords: Mathematics Subject Classification. 32V20, 35H20, 58J50 CR Sphere, Sub-Laplacian, eigenvalue

We prove that the first positive eigenvalue, normalized by the volume, of the sub-Laplacian associated with a strictly pseudoconvex pseudo-Hermitian structure θ on the CR sphere S 2n+1 ⊂ C n+1 , achieves its maximum when θ is the standard contact form.

1 n!2 n θ ∧ (dθ) n (see the next section for precise definitions). The main result of this paper is 2010

Introduction and statement of the main result

According to a classical result of Hersch [START_REF] Hersch | Quatre propriétés isopérimétriques de membranes sphériques homogènes[END_REF], given any Riemannian metric g on the 2dimensional sphere S 2 , the first positive eigenvalue λ 1 (g) of the Laplace-Beltrami operator ∆ g satisfies the estimate λ 1 (g)A(g) ≤ λ 1 (g 0 )A(g 0 )

where g 0 is the standard metric of S 2 and A(g) is the area of S 2 with respect to g. Moreover, the equality holds in [START_REF] Aribi | On the continuity of the eigenvalues of a sublaplacian[END_REF] if and only if g is isometric to g 0 . This result was extended to higher dimensional spheres by Ilias and the second author as follows (see [12, proposition 3.1]) : If a Riemannian metric g on the n-dimensional sphere S n is conformal to the standard metric g 0 , then λ 1 (g)V (g)

2 n ≤ λ 1 (g 0 )V (g 0 ) 2 n (2)
where V (g) denotes the Riemannian volume of the sphere with respect to g. Again, the equality holds in [START_REF] Aribi | Eigenvalues of the sub-Laplacian and deformations of contact structures on a compact CR manifold[END_REF] if and only if g is isometric to g 0 . The aim of the present paper is to establish a version of the estimate (2) for the first positive eigenvalue of the sub-Laplacian on the CR sphere S 2n+1 ⊂ C n+1 . Indeed, let

θ 0 = i 2 n+1 j=1
ζ j d ζjζj dζ j be the standard contact form on S 2n+1 whose kernel coincides with the Levi distribution H(S 2n+1 ) = T S 2n+1 ∩ JT S 2n+1 , where J is the complex structure of C n+1 . The set

P + (S 2n+1 ) = {f θ 0 ; f ∈ C ∞ (S 2n+1
) and f > 0} contains all pseudo-Hermitian structures on S 2n+1 whose Levi form is positive definite. Given a pseudo-Hermitian structure θ ∈ P + (S 2n+1 ), we denote by λ 1 (θ) the first positive eigenvalue of the corresponding sub-Laplacian ∆ θ , and by V (θ) the volume of S 2n+1 with respect to the volume form ψ θ = Theorem 1.1. For every pseudo-Hermitian structure θ ∈ P + (S 2n+1 ) we have

λ 1 (θ)V (θ) 1 n+1 ≤ λ 1 (θ 0 )V (θ 0 ) 1 n+1 . (3)
The equality holds in (3) if and only if there exists a CR-automorphism γ of S 2n+1 such that θ = c γ * θ 0 for some constant c > 0, or if and only if there exist p ∈ S 2n+1 and t ≥ 0 such that θ = c |cosh t + sinh t (ζ, p)| 2 θ 0 , where ( , ) denotes the standard Hermitian product of C n+1 . This result can be seen as a contribution to the program aiming to recovering the main results of spectral geometry, established for the eigenvalues of the Laplace-Beltrami operator on a compact Riemannian manifold, in the realm of CR and pseudo-Hermitian geometry. This program has motivated a lot of research in recent years and we can find significant contributions in [START_REF] Aribi | On the continuity of the eigenvalues of a sublaplacian[END_REF][START_REF] Aribi | Eigenvalues of the sub-Laplacian and deformations of contact structures on a compact CR manifold[END_REF][START_REF] Aribi | A lower bound on the spectrum of the sublaplacian[END_REF][START_REF] Aribi | Inequalities and bounds for the eigenvalues of the sub-Laplacian on a strictly pseudoconvex CR manifold[END_REF][START_REF] Barletta | On the spectrum of a strictly pseudoconvex CR manifold[END_REF][START_REF] Barletta | The Lichnerowicz theorem on CR manifolds[END_REF][START_REF] Barletta | Sublaplacians on CR manifolds[END_REF][START_REF] Chang | A zeta function associated to the sub-Laplacian on the unit sphere in C N[END_REF][START_REF] Chang | On the CR analogue of Obata's theorem in a pseudohermitian 3-manifold[END_REF][START_REF] Michael | Spectral multipliers for the Kohn sublaplacian on the sphere in C n[END_REF][START_REF] Greenleaf | The first eigenvalue of a sub-Laplacian on a pseudo-Hermitian manifold[END_REF][START_REF] Hassannezahed | Sub-laplacian eigenvalue bounds on sub-riemannian manifolds[END_REF][START_REF] Ivanov | An Obata-type theorem on a three-dimensional CR manifold[END_REF][START_REF] Ivanov | An Obata type result for the first eigenvalue of the sub-Laplacian on a CR manifold with a divergence-free torsion[END_REF][START_REF] Kokarev | Sub-Laplacian eigenvalue bounds on CR manifolds[END_REF][START_REF] Li | The sharp lower bound for the first positive eigenvalue of a sub-Laplacian on a pseudo-Hermitian manifold[END_REF][START_REF] Niu | Payne-Polya-Weinberger type inequalities for eigenvalues of nonelliptic operators[END_REF][START_REF] Raphaël | Heisenberg calculus and spectral theory of hypoelliptic operators on Heisenberg manifolds[END_REF][START_REF] Stanton | Spectral invariants of CR manifolds[END_REF].

Proof of Theorem 1.1

Let M be a connected differentiable manifold of dimension 2n + 1 ≥ 3. A CR structure on M is a couple (H(M), J) where H(M) is a 2n-dimensional subbundle of the tangent bundle T M, the so-called Levi distribution, endowed with a pseudo-complex operator J satisfying the following integrability condition :

∀X, Y ∈ Γ(H(M)), [X, Y ] -[JX, JY ] ∈ Γ(H(M)) and J ([X, Y ] -[JX, JY ]) = [JX, Y ] + [X, JY ].
Real hypersurfaces of C n+1 are the most natural examples of CR manifolds. Indeed, if M ⊂ C n+1 is such a hypersurface, then H(M) := T M ∩ J(T M) endowed with the restriction of the standard complex structure J of C n+1 , is a CR structure on M.

If (M, H(M), J) is an orientable CR manifold, then there exists a nontrivial 1-form θ ∈ Γ(T * M) such that Kerθ = H(M). Such a 1-form, called pseudo-Hermitian structure, is of course not unique. The set of pseudo-Hermitian structures consists in all the forms f θ, where f is a smooth nowhere zero function on M. To each pseudo-Hermitian structure θ we associate its Levi form L θ defined on H(M) by

L θ (X, Y ) = -dθ(JX, Y ) = θ([JX, Y ]).
The integrability of J implies that L θ is symmetric and J-invariant. The CR manifold M is called strictly pseudoconvex if L θ is definite. Of course, this condition does not depend on the choice of θ (since L f θ = f L θ ). In the sequel, we denote by P + (M) the set of all pseudo-Hermitian structures with positive definite Levi form. Every θ ∈ P + (M) is in fact a contact form which induces on M the following volume form

ψ θ = 1 2 n n! θ ∧ (dθ) n .
The associated divergence div θ is defined, for every smooth vector field Z on M, by

L Z ψ θ = div θ (Z) ψ θ .
The sub-Laplacian ∆ θ is then defined for all u ∈ C ∞ (M), by

∆ θ u = div θ (∇ H u)
where

∇ H u ∈ Γ(H(M)) is the horizontal vectorfield such that, ∀X ∈ H(M), du(X) = L θ (∇ H u, X).
The following integration by parts formula holds for any u, v ∈ C ∞ 0 (M):

ˆM (∆ θ u) v ψ θ = -ˆM L θ (∇ H u, ∇ H v) ψ θ .
Given θ ∈ P + (M), there is a unique vectorfield ξ, often called Reeb vectorfield, that satisfies θ(ξ) = 1 and ξ⌋dθ = 0. The Levi form L θ extends to a Riemannian metric on M (the Webster metric) given by

g θ (X, Y ) = L θ (X H , Y H ) + θ(X)θ(Y )
with X H = X -θ(X)ξ. The corresponding Laplace-Beltrami operator ∆ g θ is related to ∆ θ by the following formula (see [START_REF] Greenleaf | The first eigenvalue of a sub-Laplacian on a pseudo-Hermitian manifold[END_REF])

∆ θ = ∆ g θ -ξ 2 . ( 4
)
The sub-Laplacian ∆ θ is a sub-elliptic operator of order 1/2. When M is compact, it admits a self-adjoint extension to an unbounded operator of L 2 (M, ψ θ ) whose resolvent is compact (see for instance [2, Lemma 2.2]). Hence, the spectrum of -∆ θ is discrete and consists of a sequence of nonnegative eigenvalues of finite multiplicity {λ k (θ)} k≥0 with λ 0 (θ) = 0. The min-max variational principle gives

λ 1 (θ) = inf ´M u ψ θ =0 ´M |∇ H u| 2 θ ψ θ ´M u 2 ψ θ (5)
where

|∇ H u| 2 θ = L θ (∇ H u, ∇ H u).
2.1. The CR Sphere. Let S 2n+1 be the unit Sphere in C n+1

S 2n+1 = ζ = (ζ 1 , ..., ζ n+1 ) ∈ C n+1 ; j≤n+1 |ζ j | 2 = 1
endowed with its standard CR-structure. The restriction to S 2n+1 of the contact form

θ 0 = - i 2 n+1 j=1 ζj dζ j -ζ j d ζj is a pseudohermitian structure whose Reeb field ξ = i n+1 j=1 ζ j ∂ ∂ζ j -ζj ∂ ∂ ζj
generates the natural action of S 1 on S 2n+1 . Since dθ 0 is the standard Kähler form of C n+1 , the induced Levi form on H(S 2n+1 ) coincides with the standard metric of the sphere. If V p,q is the space of harmonic polynomials of bi-degree (p, q) in C n+1 , then, ξ acts on V p,q as the multiplication by i(p -q) and it can be deduced, using (4), that V p,q is an eigenspace of ∆ θ 0 on S 2n+1 with eigenvalue 2n(p + q) + 4pq (see [START_REF] Michael | Spectral multipliers for the Kohn sublaplacian on the sphere in C n[END_REF]Theorem 4.1] and [START_REF] Stanton | Spectral invariants of CR manifolds[END_REF]Proposition 4.4] for details). Therefore,

λ 1 (θ 0 ) = 2n. (6)
2.2. One-parameter groups of CR-automorphisms of the sphere. A differentiable map ϕ : M → M between two CR manifolds is a CR map if for any x ∈ M,

d x ϕ(H x (M)) ⊂ H ϕ(x) ( M) and d x ϕ • J M x = J M ϕ(x) • d x ϕ. (7)
A CR-automorphism of a CR manifold M is a diffeomorphism of M which is a CR map.

Let e n+1 = (0, • • • , 0, 1) ∈ S 2n+1 . The punctured sphere S 2n+1 \ {e n+1 } can be identified with the boundary of the so-called Siegel domain

Ω n+1 = {(z, w) ∈ C n × C : Im w > |z| 2 } ⊂ C n+1 through the CR diffeomorphism Φ : S 2n+1 \ {e n+1 } → ∂Ω n+1 given by Φ(ζ) = 1 1 -ζ n+1 (ζ 1 , • • • , ζ n , i(1 + ζ n+1 )) with Φ -1 (z, w) = 1 w + i (2iz 1 , • • • , 2iz n , w -i)
For every t ∈ R, the "dilation"

H t : ∂Ω n+1 → ∂Ω n+1 (z, w) → (e t z, e 2t w)
is a CR-automorphism of ∂Ω n+1 . We define γ t : S 2n+1 → S 2n+1 by γ t (e n+1 ) = e n+1 and, ∀ζ ∈ S 2n+1 \ {e n+1 },

γ t (ζ) = Φ -1 • H t • Φ(ζ) = 1 cosh t + sinh t ζ n+1 (ζ 1 , • • • , ζ n , sinh t + cosh t ζ n+1 ) or γ t (ζ) = 1 cosh t + sinh t ζ n+1 (ζ + (sinh t + (cosh t -1) ζ n+1 ) e n+1 ) .
Lemma 2.1. For every t, the map γ t is a CR-automorphism of S 2n+1 which satisfies

(γ t ) * θ 0 = 1 |cosh t + sinh t ζ n+1 | 2 θ 0 . Proof. Let f j (ζ) = ζ j cosh t+sinh t ζ n+1 , j ≤ n, and f n+1 (ζ) = sinh t+cosh t ζ n+1 cosh t+sinh t ζ n+1 . Then df j = dζ j cosh t + sinh t ζ n+1 - sinh t ζ j dζ n+1 (cosh t + sinh t ζ n+1 ) 2 and df n+1 = dζ n+1 (cosh t + sinh t ζ n+1 ) 2 .
Therefore

f j d fj = ζ j d ζj |cosh t + sinh t ζ n+1 | 2 - |ζ j | 2 sinh t d ζn+1 |cosh t + sinh t ζ n+1 | 2 (cosh t + sinh t ζn+1 ) , fj df j = ζj dζ j |cosh t + sinh t ζ n+1 | 2 - |ζ j | 2 sinh t dζ n+1 |cosh t + sinh t ζ n+1 | 2 (cosh t + sinh t ζ n+1 ) which gives with n j=1 |ζ j | 2 = 1 -|ζ n+1 | 2 , n j=1 f j d fj -fj df j = 1 |cosh t + sinh t ζ n+1 | 2 n j=1 ζ j d ζj -ζj dζ j + 1 -|ζ n+1 | 2 sinh t |cosh t + sinh t ζ n+1 | 2 dζ n+1 cosh t + sinh t ζ n+1 - d ζn+1 cosh t + sinh t ζn+1
.

On the other hand,

f n+1 d fn+1 -fn+1 df n+1 = 1 |cosh t + sinh t ζ n+1 | 2 sinh t + cosh t ζ n+1 cosh t + sinh t ζn+1 d ζn+1 - (sinh t + cosh t ζn+1 ) cosh t + sinh t ζ n+1 dζ n+1
Now (with

|ζ n+1 | 2 = ζ n+1 ζn+1 ), 1 -|ζ n+1 | 2 sinh t -sinh t + cosh t ζn+1 cosh t + sinh t ζ n+1 = -ζn+1 and -1 -|ζ n+1 | 2 sinh t + (sinh t + cosh t ζ n+1 ) cosh t + sinh t ζn+1 = ζ n+1 Thus, n+1 j=1 f j d fj -fj df j = 1 |cosh t + sinh t ζ n+1 | 2 n+1 j=1 ζ j d ζj -ζj dζ j that is, (γ t ) * θ 0 = i 2 n+1 j=1 f j d fj -fj df j = 1 |cosh t + sinh t ζ n+1 | 2 θ 0 .
Let p ∈ S 2n+1 be any point of the sphere and let α p ∈ U(n+1) be such that α p (p) = e n+1 . The family γ p t = α -1 p • γ t • α p is a 1-parameter group of CR-automorphisms of the sphere S 2n+1 with

γ p t (ζ) = 1 cosh t + sinh t (ζ, p) {ζ + (sinh t + (cosh t -1) (ζ, p)) p} (8) and (γ
p t ) * θ 0 = 1 |cosh t + sinh t (ζ, p)| 2 θ 0 . (9)
2.3. Preparatory lemmas.

Lemma 2.2. Let M be a strictly pseudoconvex CR manifold of dimension 2n + 1 and let θ, θ ∈ P + (M) be two pseudo-Hermitian structures with θ = f θ, f ∈ C ∞ (M). Then

ψ θ = f n+1 ψ θ (10)
Proof. From d θ = f dθ + df ∧ θ we deduce, by induction,

(d θ) n = f n (dθ) n + α n ∧ θ where α n is a differential form of degree 2n -1. Thus, θ ∧ (d θ) n = f θ ∧ (f n (dθ) n + α n ∧ θ) = f n+1 θ ∧ (dθ) n .
Lemma 2.3. Let M be a strictly pseudoconvex CR manifold of dimension 2m + 1 and let φ : M → (S 2n+1 , θ 0 ) be a CR map. Then, for every θ ∈ P + (M),

φ * θ 0 = 1 2m 2n+2 i=1 ∇ H φ i 2 θ θ (11)
where φ 1 , . . . , φ 2n+2 are the Euclidean components of φ.

Proof. Since φ is a CR map, the 1-form φ * θ 0 vanishes on H(M) which implies that there exists f ∈ C ∞ (M) such that φ * θ 0 = f θ.

(12) Differentiating, we get φ * dθ 0 = df ∧ θ + f dθ. Hence, for every X, Y ∈ H x (M), one has φ * dθ 0 (X, Y ) = f dθ(X, Y ) and, using [START_REF] Barletta | Sublaplacians on CR manifolds[END_REF],

L θ 0 (dφ(X), dφ(X)) = dθ 0 (dφ(X), J S 2n+1 dφ(X)) = dθ 0 (dφ(X), dφ(J M X)) = φ * dθ 0 (X, J M X) = f dθ(X, J M X) = f L θ (X, X).
On the other hand, since L θ 0 coincides with the standard inner product on H φ(x) (S 2n+1 ),

L θ 0 (dφ(X), dφ(X)) = 2n+2 i=1 (dφ i (X)) 2 = 2n+2 i=1 L θ (∇ H φ i , X) 2 . Thus, f L θ (X, X) = 2n+2 i=1 L θ (∇ H φ i , X) 2 .
Taking an L θ -orthonormal basis {e 1 , . . . , e 2m } of H x (M), we get

2mf = 2m j=1 2n+2 i=1 L θ (∇ H φ i , e j ) 2 = 2n+2 i=1 L θ (∇ H φ i , ∇ H φ i ) = 2n+2 i=1 ∇ H φ i 2 θ
which implies [START_REF] Ilias | Le volume conforme et ses applications d'après Li et Yau[END_REF], thanks to [START_REF] Ilias | Immersions minimales, première valeur propre du laplacien et volume conforme[END_REF].

If φ : M → R N is a map and µ is a measure on M, we denote by ´M φ dµ the vector ´M φ 1 dµ, . . . , ´M φ N dµ ∈ R N , where φ = (φ 1 , • • • , φ N ).

Lemma 2.4. Let M be a compact manifold and let µ be a measure on M such that no open set has measure zero. If φ : M → S 2n+1 is a non constant continuous map, then there exists a pair (p, t) ∈ S 2n+1 × [0, +∞) such that ˆM γ p t • φ dµ = 0.

Proof. The proof uses standard arguments (see [START_REF] Ilias | Le volume conforme et ses applications d'après Li et Yau[END_REF][START_REF] Hersch | Quatre propriétés isopérimétriques de membranes sphériques homogènes[END_REF]). We consider the map

F : S 2n+1 × (0, +∞) → B 2n+2 ⊂ R 2n+2 (p, t) → 1 V ´M γ p t
• φ dµ where V = ´M dµ and B 2n+2 is the unit Euclidean ball. Observe that (see [START_REF] Chang | A zeta function associated to the sub-Laplacian on the unit sphere in C N[END_REF]), ∀p ∈ S 2n+1 , γ p 0 is the identity map while, ∀ζ = -p, γ p t (ζ) tends to p as t → +∞. Consequently, F (•, 0) = 1

V ´M ϕ dµ is a constant map and F (•, t) tends to the identity of S 2n+1 as t → +∞. Such a map F is necessarily onto which implies that the origin of R 2n+2 belongs to its image.

2.4. Proof Theorem 1.1. Let θ ∈ P + (S 2n+1 ) be a strictly pseudoconvex pseudo-Hermitian structure. We apply Lemma 2.4 to the identity map of S 2n+1 and the measure induced by ψ θ to obtain the existence of a pair (p, t) ∈ S 2n+1 × [0, +∞) such that ˆS2n+1 γ p t ψ θ = 0.

For simplicity, we write γ for γ p t . The Euclidean components γ 1 , . . . γ 2n+2 of γ satisfy ´S2n+1 γ j ψ θ = 0. Hence, applying the min-max principle (5), we get for every j ≤ 2n + 2,

λ 1 (θ) ˆS2n+1 γ 2 j ψ θ ≤ ˆS2n+1 ∇ H γ j 2 θ ψ θ .
Summing up we obtain, with j≤2n+2 γ 2 j = 1,

λ 1 (θ)V (θ) ≤ j≤2n+2 ˆS2n+1 ∇ H γ j 2 θ ψ θ (13)
≤ 2n

ˆS2n+1 1 2n ∇ H γ j 2 θ n+1 ψ θ 1 n+1 V (θ) 1-1 n+1 . (14)
Using Lemma 2.3, we see that, since γ is a CR map from (S 2n+1 , θ) to (S 2n+1 , θ 0 ),

γ * θ 0 = 1 2n j≤2n+2 ∇ H γ j 2 θ θ (15)
which gives, thanks to Lemma 2.2

ψ γ * θ 0 = 1 2n j≤2n+2 ∇ H γ j 2 θ n+1 ψ θ .
Thus λ 1 (θ)V (θ) ≤ 2nV (γ * θ 0 )

1 n+1 V (θ) 1-1 n+1
. Since V (γ * θ 0 ) = V (θ 0 ), we finally get λ 1 (θ)V (θ)

1 n+1 ≤ 2nV (θ 0 ) 1 n+1
which proves the inequality of the theorem thanks to [START_REF] Barletta | On the spectrum of a strictly pseudoconvex CR manifold[END_REF]. Now, if the equality holds in [START_REF] Aribi | A lower bound on the spectrum of the sublaplacian[END_REF], this means that we have equality in the Cauchy-Schwarz inequality used in [START_REF] Hassannezahed | Sub-laplacian eigenvalue bounds on sub-riemannian manifolds[END_REF]. Thus, j≤2n+2 ∇ H γ j 2 θ is constant and, thanks to [START_REF] Hersch | Quatre propriétés isopérimétriques de membranes sphériques homogènes[END_REF], θ is proportional to γ * θ 0 which takes the form given by [START_REF] Chang | On the CR analogue of Obata's theorem in a pseudohermitian 3-manifold[END_REF]. Conversely, it is clear that when γ is a CR-automorphism of the sphere then λ 1 (γ * θ 0 ) = λ 1 (θ 0 ) = 2n and V (γ * θ 0 ) = V (θ 0 ). Hence, the equality holds in (3).