Audric Drogoul

Gilles Aubert

Bruno Cessac

Pierre Kornprobst

A nonconvex variational approach for receptive field estimation

Keywords: Neurosciences, inverse problem, receptive field estimation, variational approach, nonsmooth and nonconvex optimization, Kurdyka-Lojasiewicz Neurosciences, problème inverse, estimation de champs récepteur, approche variationnelle, optimisation

Determining the receptive field of a visual sensory neuron is a first but crucial step towards the characterization of neurons response to local spatio-temporal stimuli. Existing methods are based on convex optimization methods neglecting biophysical constraints of neurons (bounded firing rate), and they are relatively poor in terms of accuracy and running time. We propose a new method to estimate receptive fields by a nonconvex variational approach, thus relaxing the simplifying and unrealistic assumption of convexity made by standard approaches. The method consists in studying a relaxed discrete energy minimized by a proximal alternating minimization algorithm. We compare our approach with the classical spike-triggered-average technique on simulated data, considering a typical retinal ganglion cell. Results show a high improvement in term of accuracy and convergence with respect to the duration of the experiment.

Introduction

This paper deals with the estimation of receptive fields of individual visual neurons (e.g., ganglion cells in the retina). Knowing the receptive field of a particular neuron allows to characterize the relation between stimulus (images) and the neuron response (action potentials, also called spikes). It is a first but crucial step to understand which region of the visual field and which stimuli a neuron is sensitive to. Thus, biologists always allow time for estimating receptive field in their experimental protocole. * A classical view is that visual sensory neurons respond to specific local spatio-temporal patterns, and that their response can be characterized by a spatio-temporal convolution kernel called the receptive field, followed by a static nonlinearity and stochastic (Poisson-like) mechanisms of spikes generation. The static nonlinearity is in general non convex (e.g., sigmoid shape) since firing rate is bounded because neurons have a refractory period. Mathematically, this description corresponds to the so-called linear-nonlinear Poissonian (LNP) model [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF][START_REF] Simoncelli | The Cognitive Neurosciences, 3rd edition, chapter Characterization of neural responses with stochastic stimuli[END_REF][START_REF] Pillow | Prediction and decoding of retinal ganglion cell responses with a probabilistic spiking model[END_REF][START_REF] Paninski | Statistical models for neural encoding, decoding, and optimal stimulus design[END_REF]. LNP models can simulate the spiking activity of ganglion cells (or cortical cells) in response to synthetic or natural images [START_REF] Carandini | Do we know what the early visual system does[END_REF], and they voluntarily ignore the detailed neuronal mechanisms. These functional model are widely used by experimentalists to characterise the cells that they record, map their receptive field and characterise their spatiotemporal feature selectivities [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF].

The LNP model is illustrated in Fig. 1: Given a visual stimulus s : (x, t) ∈ Ω × [0, T ] → R, where Ω ⊂ R 2 is the spatial domain and T > 0 is the duration of the experiment, a LNP neuron generates a sequence of n(T ) spikes times {t i } 1≤i≤n(T ) such that {t i } 1≤i≤n(T ) is generated by a Poisson process of rate r(t) = f ((s × u)(t)) , (P)

where f is a nonlinear function and u : Ω d := Ω × [-d, 0] -→ R is the so-called receptive field which corresponds to the linear part of the processing, where d > 0 is the length of its temporal support: the response of the neuron at time t depends on the history of the stimulus up to time td back in the past. The operator × is defined by

[s × u](t) = 0 -∞ Ω
s(x, t + τ )u(x, τ )dx dτ.

(1)

Assuming Hypothesis (P), the problem of estimating the receptive field u(x, t) is an inverse problem: given a stimulus s(x, t) of duration T and the n(T ) spikes {t i } 1≤i≤n(T ) , recover the unknown receptive field u(x, t). Here the stimulus is a sequence of white noise images, which is the classical stimulus used for the spike-triggered averaged approach (STA, see [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF]).

This inverse problem can be formulated using a Bayesian approach formulated in a discrete setting. Considering a discretization of the receptive field u N ∈ R N , we search u N as maximum likelihood or equivalently (see, e.g., [START_REF] Paninski | Maximum likelihood estimation of cascade point-process neural encoding models[END_REF][START_REF] Pillow | Likelihood-based approaches to modeling the neural code. Bayesian brain: Probabilistic approaches to neural coding[END_REF][START_REF] Park | Receptive Field Inference with Localized Priors[END_REF]):

inf u N -log ρ(u N |{t i } 1≤i≤n(T ) ) = inf u N -log(ρ({t i } 1≤i≤n(T ) |u N )) -log(ρ(u N )) , (2)
where ρ(u N |{t i } 1≤i≤n(T ) ) is the density probability of the random variable u N given {t i } 1≤i≤n(T ) , ρ({t i } 1≤i≤n(T ) |u N ) is the density probability of the random variable {t i } 1≤i≤n(T ) given u N , and ρ(u N ) is the density probability of the random variable u N . Then, assuming that ρ(u N ) is a distribution of the form ρ(u N ) = e -J(u N ) , and passing to the limit when N → ∞ in [START_REF] Beck | Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems[END_REF] we get that the continuous problem is (see Appendix, Sect. A):

inf u E(u) = ψ(s × u) + J(u), (3)
where ψ(s × u) is the data fidelity term with ψ defined by

ψ(z) = T 0 f (z(τ ))dτ - n(T ) i=1 log(f (z(t i ))), (4)
and J(u) is the prior term to infer qualitative properties to the solution, according to what is known about the general shape of a receptive field.

To solve [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF], it is classical to assume that the nonlinearity f is a ramp function or an exponential. In that case, one has to solve a convex problem which can be done by classical methods [START_REF] Paninski | Maximum likelihood estimation of cascade point-process neural encoding models[END_REF][START_REF] Pillow | Likelihood-based approaches to modeling the neural code. Bayesian brain: Probabilistic approaches to neural coding[END_REF][START_REF] Mcfarland | Inferring Nonlinear Neuronal Computation Based on Physiologically Plausible Inputs[END_REF]. On the opposite, in this paper we consider the case when the nonlinearity f can be a sigmoid, which is more realistic from a physiological point of view: it models the fact that the neuron has a bounded firing rate (see [START_REF] Gerstner | Spiking Neuron Models[END_REF], Sect.5.2.3). However, since f is nonconvex, the data fidelity term ψ(s × u) also becomes nonconvex.

The prior term J(u) is chosen according to the properties we want to impose on the solution. The first property is that u should be localized in space and time since neurons are sensitive to a particular region of the visual field. This can be imposed by a sparsity constraint term. Here we choose to use a convex relaxation of the sparsity and we penalise the L 1 -norm of u [START_REF] Candes | Enhancing sparsity by reweighted 1 minimization[END_REF]. The second property is that u should be smooth. We propose that u belongs to the space BV 2 that contains piecewise linear functions. Note that in [START_REF] Bergounioux | A Second-Order Model for Image Denoising[END_REF], this regularity constraint has been used for image restoration. It allows to recover functions with fast smooth variations more precisely than using a simple ∇u 2 L 2 -regularity constraint. The space BV 2 , is defined by

BV 2 (O) = u ∈ W 1,1 (O), ∂u ∂x i ∈ BV (O) for i ∈ {1, ..., N } ,
where O ⊂ R N with N ∈ N is a regular domain, W 1,1 (O) and BV (O) stands respectively for the Sobolev space and the space of functions of bounded variations. On BV 2 , the following norm can be defined:

|u| BV2(O) ≡ O |D 2 u| = sup O < ∇u, div(ϕ) >, ϕ ∈ C 1 0 (O) N ×N , ϕ L ∞ (O) ≤ 1 ,
being the total variation of the gradient where C 1 0 (O) is the space of functions C 1 (O) with compact support in O [START_REF] Bergounioux | A Second-Order Model for Image Denoising[END_REF][START_REF] Demengel | Fonctions à hessien borné[END_REF]. So the prior term will be defined by

J(u) = λ u L 1 (Ω d ) + µ|u| BV2(Ω d ) , (5)
where λ, µ > 0 are weights associated respectively to the sparsity and to the regularity of u. This paper is concerned with the study of problem (3) with ψ and J defined respectively by (4) and [START_REF] Berdondini | Active pixel sensor array for high spatio-temporal resolution electrophysiological recordings from single cell to large scale neuronal networks[END_REF]. In Sect. 2, we study the discrete version of problem [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF]. We introduce a relaxed problem which allows to compute an approximation of the solution. We make a theoretical study of the relaxed problem and propose and alternated minimizing algorithm converging toward a critical point of the relaxed energy. In Sect. 3, we test the approach on simulated data to provide a quantitative evaluation with comparisons to the classical spike triggered averaged technique (STA, see [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF]) and we compare our approach with the classical variational approach using as non linearity a convex function close to the nonlinear sigmoid function describing the firing rate of a real neuron. Section 4 is a discussion and conclusion section. Let us also mention several appendices presenting more technical developments so that the paper be self-contained.

The discrete problem: well-posedness and relaxation

In this section we study the discrete problem. Section 2.1 gives some mathematical preliminaries. Section 2.3 is about the well-posedness. Section 2.4 presents the main algorithm which is an alternating minimization procedure followed by sections describing algorithms to compute solutions of each sub-problem.

Preliminaries

We first recall some classical concepts of variational analysis which will be useful in the sequel. Let φ : R n -→ R∪{+∞} a proper lower semicontinuous function. We denote as usual dom(φ

) = {x ∈ R n , φ(x) < +∞}, epi(φ) = {(x, λ) ∈ R n , φ(x) ≥ λ} and graph(φ) = {(x, λ) ∈ R n , φ(x) = λ}
respectively the domain, the epigraph and the graph of φ.

Notions of subdifferentials

We recall that the Fréchet subdifferential of φ at x ∈ dom(f ) is defined by (see [START_REF] Rockafellar | Variational Analysis, volume 317 of Grundlehren der mathematischen Wissenschaften[END_REF], Definition 8.3)

∂ F φ(x) = x ∈ R n , lim inf y→x y =x 1 x -y (φ(y) -φ(x) -x , y -x ) ≥ 0 .
The limiting subdifferential of φ at x ∈ dom(φ) is denoted and defined by

∂φ(x) = {x ∈ R n ∃x n → x, φ(x n ) → φ(x), x n ∈ ∂ F φ(x n ) → x }.
Let us remark that

∂ F φ ⊂ ∂φ,
where the first set is convex and closed while the second is closed (see [START_REF] Rockafellar | Variational Analysis, volume 317 of Grundlehren der mathematischen Wissenschaften[END_REF], Theorem 8.6).

We say that φ is convex if its epigraph is convex, and semiconvex if there exists a constant c ≥ 0, such that φ(z)

+ c 2 z 2 is convex. If φ is convex, for x ∈ dom(φ) we have ∂ F φ(x) = ∂φ(x) = ∂φ(x) def := {x ∈ R n , φ(x) -φ(x ) ≥ x , x -x ∀x ∈ R n }.
These equalities remains true for semiconvex functions (see [START_REF] Rockafellar | Variational Analysis, volume 317 of Grundlehren der mathematischen Wissenschaften[END_REF], Example 8.8).

We say that x is a critical point of φ if and only 0 ∈ ∂φ(x).

Finally, we say that φ : R n → R is separable if it can be written as φ(x) = n i=1 φ i (x i ) where φ i are real functions.

Inria

Notion of proximal operator

The following mapping introduced by Moreau [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] is very useful in non-smooth optimisation, and is a generalization of the notion of projector. The proximal operator of φ is denoted and defined by:

prox φ (x) = arg min y∈R n φ(y) + 1 2 x -y 2 2 ∀x ∈ R n , (6)
and it can be shown that

prox φ = I + ∂φ -1
.

Remarks 2.1. Some properties of the proximal operator. 2 is convex for all x, then prox φ is univalued. (ii) If φ is not convex, the proximal operator prox φ can be multivalued.

(i) If φ + 1 2 (. -x)
(iii) Let λ > 0, the function φ λ : x → inf y {φ(y) + 1 2λ x -y 2 2
} is called the Moreau envelope or Moreau-Yoshida regularization of φ. The Moreau envelope can be considered as a smooth approximation of φ. Note that φ λ (x) ≤ φ(x)∀x.

(iv) Let us consider a shifted and scaled proper function φ(x) = φ(ax + b), then its proximal operator is given by prox φ(x) = (prox a 2 φ (ax + b)b)/a.

Notion of conjugate function

The Legendre-Fenchel transform of φ is defined by (see [START_REF] Rockafellar | Variational Analysis, volume 317 of Grundlehren der mathematischen Wissenschaften[END_REF], p. 473)

φ (x) = sup y x, y -φ(y) ∀x ∈ R n . ( 7
)
Function φ is called the conjugate of function φ, and φ = (φ ) is the biconjugate of function φ. For a subset K ⊂ R n we define the indicator function i K and the support function δ K by

i K (x) = 0 if x ∈ K, +∞ otherwise, and δ K (x) = sup x∈K x, y . (8)
The link between δ K and i K is that they are mutually conjugate when K is convex and closed. Let us recall that all functions φ lower semicontinous, convex and one-homogeneous (i.e., φ(λx) = λφ(x) for λ > 0) verify

φ (x) = δ K φ (x) with K φ = {x ∈ R n , ∀y ∈ R n y, x ≤ φ(y)}. (9)

Problem definition

In this section, we study a discrete version of (3). The stimulus is a sequence of N t images, each one presented during a period of ∆t so that the duration of the experiment is

T = N t ∆t. Each image is of size N x × N y pixels. Receptive field is of size N x × N y × D
where D is a fixed depth in time. We introduce the real vector spaces X , Y and Z so that:

u ∈ X = R Nx×Ny×D , s ∈ Y = R Nx×Ny×Nt , z ∈ Z = R Nt , (10)
endowed with the scalar product ., . and the associated norm . 2 = ., . . We denote by . 1 the l 1 -norm.

Given these notations, if ξ = (ξ i ) 1≤i≤Nt is the number of spikes per time step (in [i∆t, (i + 1)∆t[) and if we assume that ∆t = 1 without loss of generality (up to change f by ∆tf ), then the data fidelity term (4) can be rewritten as

ψ ξ (z) = Nt i=1 f (z i ) - n(T ) i=1 log(f (z ki )),
where k i is the index of the small interval containing the spiking time t i , so that

ψ ξ (z) = Nt i=1 ψ ξi (z i ), with ψ ξi (z i ) = f (z i ) -ξ i log(f (z i )). (11)
The function ψ ξi (z) has a limited domain of definition denoted by Z + ξi equal to R if ξ i = 0 (no spikes), and {f > 0} if ξ i > 0.

We impose the following hypotheses on the sigmoid function f :

Hypothesis 1. f is C 0 (R)
and is bounded from below.

Hypothesis 2. (i) f is definable in the log-exp structure (see Appendix C for definition and properties around these notions).

(ii) f is semi-convex, i.e., there exists ρ > 0 such that

f (z) + ρ 2 z 2 is convex. (iii) z → -log(f (z)) is convex on the set {f > 0}.
Remark 2.1. Hypothesis 2(i) is required by Theorem 2.1 concerning the convergence of the algorithm (see Sect. 2.4). Note that the assumption to be definable in the log-exp structure is not very restrictive. For example, the nonlinearities generally chosen are C 0 (R) and analytic on R except at a finite number of points (see, e.g., [START_REF] Paninski | Maximum likelihood estimation of cascade point-process neural encoding models[END_REF]), so that they verify this hypothesis. Concerning Hypothesis 2(ii), remark that it is verified by convex functions, so that all results coming next will hold in the convex case too.

In practice, we will choose f as a non decreasing piecewise cubic sigmoid function varying on [θ 1 , θ 2 ], equals to zeros for z < θ 1 and equals to a constant c > 0 on [θ 2 , +∞[. Figure 2 shows the shape of the data fidelity term ψ ξ (.) for this kind of sigmoid. One can easily verify that the data fidelity term is nonconvex when 0 ≤ ξ < c.

✓ 1 ✓ 2 0 c ✓ 1 ✓ 2 f 1 (⇠) c ⇠log(c) ✓ 1 ✓ 2 c ⇠log(c) c ⇠log(c) c ⇠log(c) ⇠ (.) ⇠ (.) ⇠ (.) ⇠ = 0 0 < ⇠ < c ⇠ c 0 0 0

Inria

The prior term is defined by

J(u) = λ u 1 + µ Hu 1 , (12)
where H : X → X ν is the Hessian matrix operator (ν = 9 in 3D).

To discretize the Hessian operator, we use symmetric boundary conditions: (Hu) i,j,k = ((Hu) i,j,k )

p,q for 1 ≤ p, q ≤ 3 with

(Hu) 1,1 ijk =      u i+1,j,k -2u i,j,k + u i-1,j,k if 1 < i < N x , u i+1,j,k -u i,j,k if i = 1, -(u i,j,k -u i-1,j,k ) if i = N x , (Hu) 1,2 ijk =      u i,j+1,k -u i-1,j+1,k -u i,j,k + u i-1,j,k if 1 < i < N x and 1 < j < N y , 0 if i = 1, 0 if j = N y ,
Other derivatives can be obtained by permutation of the indices.

Denoting by × d the discretization of × we get the discrete problem associated to (3):

inf u∈X E(u) = ψ ξ (s × d u) + i Z + ξ (s × d u) + λ u 1 + µ Hu 1 , (13)
where i Z + ξ denotes the characteristic function associated to

Z + ξ = Nt i=1 Z + ξi .
In [START_REF] Carandini | Do we know what the early visual system does[END_REF], the data fidelity term only depends on s × d u and when minimizing it w.r.t. u, it is difficult to guarantee that s × d u remains in the domain of definition of ψ ξ . To overcome this problem, we propose to introduce an auxiliary variable z ∈ Z and solve the following relaxed formulation:

inf z∈Z,u∈X E α (z, u) = ψ ξ (z) + i Z + ξ (z) + α 2 s × d u -z 2 2 + λ u 1 + µ Hu 1 , (14)
where the term in α penalizes the difference between z and s × d u. Another interest of this relaxed problem is that now the problem in z containing the nonconvexity is separable (it leads to N t one-dimensional independent problems). This relaxation is close to augmented Lagrangian for alternating direction method of multipliers (ADMM) algorithms [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF] used in a convex context (this later just introduces another dual variable to guarantee the equality constraint at critical point and consists in two descents and one ascent of gradient). Let us notice that [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF] gives an approximation of solutions of (13) when α → ∞ while ADMM algorithms with augmented Lagrangian give exact solutions for a convex energy E. In this paper we consider the general case where E α is not convex and we will apply existing generic algorithms to compute critical points of E α which contain solutions of (14).

Well-posedness

This section is about existence and uniqueness of solutions of problems ( 13) and ( 14).

Proposition 2.1 (existence). Assuming Hypothesis 1, then 1. Let λ > 0, the energy E(u) (13) admits at least one minimizer.

2. Let α > 0 and λ > 0, the energy E α (z, u) (14) admits at least one minimizer.

Proof. We only prove that E α (z, u) admits a minimizer for α > 0 and λ > 0, the existence of minimizer of E(u) being similar. Thanks to compactness of bounded sequences in R n and lower semicontinuity of E α (z, u) (E α is even continuous), the existence question reduces to the sufficient question of coerciveness. Let (z n , u n ) n , a sequence such that E α (z n , u n ) be bounded by a constant M > 0. As ψ ξ (z) is bounded from below, we deduce that u n 1 and z ns × d u n are bounded. From that it is clear that z n is also bounded. Thus by extracting subsequences and using the lower semicontinuity of E α , we easily deduce that there exists (z, ū) such that

E α (z, ū) = inf E α . Proposition
z i → f (z i ) + α 2 (z i -x) 2 is strictly convex for each x ∈ R if α > ρ.
From the separability of E α (z, u) w.r.t z and Hypothesis 2(iii) we deduce that z → E α (z, u) is strictly convex. The first term z → Ψ ξ (z) is convex, the second term

(z, u) → s × d u -z 2
2 also, and the third term u → J(u) as well. Hence, we get that E α (z, u) is convex and we deduce that critical points are global minima. Indeed, we have

E α (z, u) -E α (z, ū) ≥ ∂E α (z, ū).(z -z, u -ū) = 0, ∀ (z, u) ∈ Z × X ,
for each critical points (z, ū). Remark 2.2. Even if we assume that f is strictly convex, the energy

E α (z, u) is in general not strictly convex. The reason is that u → α 2 z -s × d u 2 2 + φ(Ku) with K := I × H : X → X × X ν and φ(v) = i |v 1i | + |v 2i |
is in general not strictly convex since φ is not strictly convex and u → s × d u is in general not injective.

The following proposition establishes the link between solutions of the relaxed problem and solutions of the initial one.

Proposition 2.3. Let α be a positive integer and let {(z α , u α )} α≥1 be a family of solution of problem [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF] then E α (z, u) is bounded independently of α and (i) All its cluster points are couples (s× d ū, ū) such that ū is a solution of (13).

(ii) The infimum converges:

inf (z,u) E α (z, u) = E α (z α , u α ) -→ α→+∞ inf u E(u). (iii) If E(u) admits a unique minimizer ū, then (z α , u α ) -→ α→+∞ (s × d ū, ū).
Proof. Let be ū a minimizer of the initial energy E(u) [START_REF] Carandini | Do we know what the early visual system does[END_REF] and (z α , u α ) a minimizer of the relaxed energy E α (z, u) [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF]. Let α 0 > 0, and α ≥ α 0 , we have that

E α0 (z α , u α ) ≤ E α (z α , u α ) ≤ E α (s× d ū, ū) = E(ū). (15)
As in the proof of Prop. 2.1, we deduce that {(z α , u α )} α>0 is bounded independently of α and that it converges up to a subsequence to (z, ũ) ∈ Z + ξ × X . From the definition of E α (z, u) and since ψ ξ (z) and J(u) are bounded from below we deduce that z α -s× d u α -→ α→+∞ 0 and that z = s× d ũ. Now let us prove that ũ is a minimizer of E(u). It is easily seen that the quantity m α = E α (z α , u α ) is non-decreasing and bounded from above by E(ū) so it converges to its supremum. By uniqueness of the limit we have that lim α m α = E(ũ) and by passing to the limit into [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF], we get that E(ũ) = E(ū). Hence ũ is a minimiser of E(u).

Inria

Main algorithm

To compute a solution of ( 14), we propose the proximal alternating minimization introduced in [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF] in a more general context. Given an initial condition u 0 ∈ X , the algorithm consists in the following two steps:

z (k+1) ∈ arg min z E α (z, u (k) ) + 1 2β (k) z -z (k) 2 2 , (16a)
u (k+1) ∈ arg min u E α (z (k+1) , u) + 1 2γ (k) u -u (k) 2 2 , (16b)
where β (k) , γ (k) are sequences of parameters belonging to [r -, r + ] with 0 < r -< r + for all k ≥ 0 (note that this is the only condition on these parameters to obtain convergence). The quadratic terms in [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] are necessary to show the convergence and it can be shown that one minimization step leads to a decrease of the energy such that

E α (z (k+1) , u (k+1) ) + 1 2β (k) z (k+1) -z (k) 2 2 + 1 2γ (k) u (k+1) -u (k) 2 2 ≤ E α (z (k) , u (k) ),
which gives an interpretation of the qualitative role of parameters β (k) , γ (k) for convergence. Thus, the algorithm ( 16) consists at each iteration to compute the proximal mapping (up to a multiplicative constant) of z → E α (z, u (k) ) and of u → E α (z (k+1) , u) successively. From a numerical point of view, this algorithm is easy to handle:

• If f is a piecewise cubic function verifying Hypotheses 1 and 2, when ξ i = 0, then problem (16a) can be solved analytically thanks to its separability: it is equivalent to compute the proximal operator of f up to a multiplicative constant depending on α and β (k) . When α + 1 β (k) > ρ the problem is strictly convex so that the proximal operator is univalued. Otherwise, in the general case set by Hypotheses 1 and 2, we use a Newton algorithm.

• Problem (16b) is strictly convex. It can be solved by regularising the non differentiable and non separable term Hu 1 and by using a standard forward-backward proximal algorithm of Nesterov type [START_REF] Beck | Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems[END_REF]. Another possibility is to rewrite the non differentiable term of (16b) by using the Fenchel transform [START_REF] Bolte | The Łojasiewicz Inequality for Nonsmooth Subanalytic Functions with Applications to Subgradient Dynamical Systems[END_REF] and to introduce a dual variable, inspired by [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF].

Theorem 2.1. (convergence) Assume Hypotheses 1 and

2(i), Algorithm (16a)-(16b) generates a sequence (z (k) , u (k) ) k which converges to a critical point (z, ū) of E α (z, u) (i.e., 0 ∈ ∂E α (z, ū)) and the sequence (z (k) , u (k) ) verifies ∞ k=1 z (k) -z (k-1) + u (k) -u (k-1) < ∞.
Proof. To prove convergence, we apply Theorem 9 from [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF]. The difficult point is to check that E α verifies the Kurdyka-Lojasiewicz (KL) property which is a notion introduced in 1963 by Lojasiewicz for analytical functions and extended by Kurdyka in 1998 to functions definable in an o-minimal structure. This is detailed in Appendix C.

Remark 2.3. The symbol ∈ can be replaced by = in (16b), since the subproblem is strictly convex for α > 0. If we furthermore assume Hypothesis 2(ii)-(iii) then we can do as well in (16a). What is remarkable is that the sequence (z (k) , u (k) ) converges to a critical point of E α even if each subproblem of (16) has a non unique solution, and from our knowledge this theorem is the only one that establishes the convergence of the whole sequence to a critical point, just assuming that the energy verifies the KL property (Definition C.1) and some classical assumptions to guarantee existence of minimizer(s).

Corollary 2.1. Assume that f verifies Hypothesis 1 and 2 and moreover that f is convex, then the sequence (z (k) , u (k) ) k converges to a global minimum of E α .

Proof. Under these assumptions and thanks to Prop. 2.2, we get that E α is convex so critical points are global minima which ends the proofs by applying Th. 2.1.

2.5 Solving problem (16a)

Using proximal operator definition

We remark that E α (z, u) is separable w.r.t z. Hence minimizing it w.r.t z at fixed u is equivalent to minimize N t one dimensional energies:

z (k+1) i = arg min zi∈R ψ ξi (z i ) + i Z + ξ i (z i ) + β (k) α + 1 2β (k) z i - αβ (k) (s× d u (k) ) i + z (k) i β (k) α + 1 2 ,
where ξ i is the number of spikes in the i-th time interval and ψ ξi defined in [START_REF] Butts | Temporal Precision in the Visual Pathway through the Interplay of Excitation and Stimulus-Driven Suppression[END_REF] is the data fidelity term. We rewrite z (k+1) i by using the proximal operator ( 6) associated to the i-th component of the data fidelity term:

z (k+1) i = prox h (k) ψ ξ i αβ (k) (s × d u (k) ) i + z (k) i β (k) α + 1 , (17)
where

h (k) = β (k) β (k) α+1
and where we shorten the notation

ψ ξi + i Z + ξ i
by ψ ξi . The function ψ ξi depends on the number of spikes ξ i and is defined from the sigmoid function. Generally the sigmoid function is approximated by an exponential based function (see Hypothesis 2(i) and Fig. 3(a)). The problem is that finding an analytical expression of z

(k+1) i
in that case is a hard task. As a consequence, we propose two approximations for which analytical expressions can be derived. The first approximation is defined by a piecewise linear approximation (see Fig. 3(b) and Appendix D) but in this case problem (16a) is nonconvex independently of α and β (k) (the proximal operator is multivalued). The second approximation is defined by a piecewise cubic sigmoid function (see Fig. 3(c) and Appendix E). It is the simplest choice verifying Hypotheses 1 and 2 which enables to compute z (k+1) i (17) analytically when there is no spike (ξ i = 0) and numerically otherwise. This fact is interesting because spikes are sparsely distributed.

In the two following sections we study [START_REF] Demengel | Fonctions à hessien borné[END_REF] for ξ i = 0 (no spike in the i-th interval) and ξ i > 0 (at least one spike).

Case

ξ i = 0 (no spike)
In this case, [START_REF] Demengel | Fonctions à hessien borné[END_REF] reduces to the computation of the proximal operator of γf with γ = h (k) . Thanks to Remark 2.1-(iv), without loss of generality, we focus on a normalized sigmoid function. Considering different nonlinear functions f , Table 1 summarizes their properties, if they verify Hypotheses 1 and/or 2 and their proximal operator. We give results for the three nonconvex functions described in Fig. 3 but also for convex functions for which the same method applies (see numerical results in Sect. 3). Details about how to estimate proximal operators are given in Appendices D and E for piecewise linear and piecewise cubic sigmoid cases respectively. [START_REF] Rockafellar | Variational Analysis, volume 317 of Grundlehren der mathematischen Wissenschaften[END_REF]. (c) Piecewise cubic approximation [START_REF] Su | A differential equation for modeling nesterov's accelerated gradient method: Theory and insights[END_REF].

Inria -1 2 0 1 2 0 1 1 2 -1 2 0 1 2 0 1 1 2 -1 2 0 1 2 0 1 1 2 (a) (b) (c)
Nonlinear function f (x) Properties Proximal operator prox γf (x)      0 if x < -1 2 1 2 + x if |x| ≤ 1 2 1 if x > 1 2 Hyp. 1 Hyp. 2(i)-(iii) nonconvex          x if x ≤ -1 2 -1 2 if -1 2 ≤ x ≤ γ x -γ if -1+2γ 2 ≤ x ≤ 1+γ 2 x if x ≥ 1 2      0 if x ≤ -1 2 1 2 + 3 2 x -2x 3 if |x| ≤ 1 2 1 if x ≥ 1 2
Hyp. 1 and 2 nonconvex

     x if x ≤ -1 2 y -(x) if |x| ≤ 1 2 x if x ≥ 1 2 y -(x)
given by ( 42)

multivalued if γ > 1 6 1 1 + e -x
Hyp. 1 and 2 nonconvex

No analytical expression

Use a Newton algorithm

0 if x ≤ -1 2 1 2 + x if x ≥ -1 2 Hyp. 1 and 2 convex      x if x ≤ -1 2 -1 2 if -1 2 ≤ x ≤ -1+2γ 2 x -2γ if x ≥ -1+2γ 2 0 if x ≤ -1 2 1 2 + 2x + 2x 2 if x ≥ -1 2
Hyp. 1 and 2 convex

x if x ≤ -1 2 x -2γ 4γ + 1 if x ≥ -1 2 log(1 + e x )
Hyp. 1 and 2 convex

No analytical expression

Use a Newton algorithm Nonlinear function f (x) Properties Proximal operator prox γψ ξ (x)

     0 if x < -1 2 1 2 + x if |x| ≤ 1 2 1 if x > 1 2 Hyp. 1 Hyp. 2(i)-(iii) nonconvex      y + (x) if x ≤ x 1 1 2 if x 1 ≤ x ≤ 1 2 x if x ≥ 1 2 y + (x)
given by ( 40)

x 1 = γ(1 -ξ) + 1 2 multivalued if ξ < 1 Non decreasing non nega- tive and C 0 , piecewise C 1 -    ∈ ] -1 2 , max(x, f (-1) (ξ)] if f (-1) (ξ) = ∅ x otherwise
Use a Newton algorithm introducing a regularizing parameter † or use a Golden search algorithm.

Table 2: Proximal operators of γψ ξ (with ψ ξ (x) = f (x)ξlog(f (x))) for different nonlinear functions f (see Sect. 2.5.3).

Case ξ i > 0 (presence of spikes)

In this case, (17) expands as

z (k+1) i = arg min y∈Z + ξ i f (y) -ξ i log(f (y)) + 1 2h (k) (y -a i,k ) 2 , ( 18
)
with a i,k =

αβ (k) (s× d u (k) )i+z (k) i β (k) α+1
and

h (k) = β (k) αβ (k) +1 .
As stated in Prop. 2.2, assuming Hypothesis 2(i)-(ii), problem [START_REF] Dries | Tame topology and o-minimal structures[END_REF] is strictly convex for α > ρ -1 β (k) , hence it admits a unique solution that is the proximal operator of h (k) ψ ξ evaluated at a i,k . If we model the nonlinearity by a piecewise linear sigmoid function (which does not satisfy Hyp. 2(ii)), we get an analytical expression of the proximal operator of γψ ξ for some γ > 0 which is given in Table 2. Otherwise, to compute efficiently an approximation of z (k+1) i , either we use a Newton algorithm by regularizing the term -log(f (y)) by -log(f (y) + η) with η > 0 small enough or we use a Golden search algorithm by using the fact that the solution belongs to ]θ 1 , z max [ with z max = max(a i,k , f -1 (ξ)) with the convention f -1 (x) = ∅ if x is not in the range of f as summed up in Table 2.

Remark 2.4. Let us remark that if we know the proximal operator associated to ψ f ξ = f -ξlog(f ) for a given non linearity f we can deduce the one associated to ψ g ξ for whatever g(x) = cf (ax + b) thanks to this relation

ψ g ξ (x) = cψ f ξ c (ax + b) -ξlog(c),
and Remark 2.1(iv). † In Newton algorithm, if f can be null, one can regularize log(f ) by changing it to log(f + η) with η > 0 small enough, so that dom(ψ ξ ) = R.

Inria

Conclusion on a good choice of sigmoid

By modelling the sigmoid function by a C 1 (R) piecewise cubic function, Hypotheses 1 and 2 are verified and problem ( 14) is convex for α large enough (Proposition 2.2). In this case we also get an analytical expression of the proximal operator associated to the data fidelity term [START_REF] Butts | Temporal Precision in the Visual Pathway through the Interplay of Excitation and Stimulus-Driven Suppression[END_REF] in the most frequent case which is the case of no spike has been emitted in the small i-th temporal bin (ξ i = 0). It remains true for some other nonlinearities enumerated in Table 1. The choice of a piecewise linear sigmoid function enables to have an analytical expression of the proximal operator associated to the data fidelity term γψ ξ [START_REF] Butts | Temporal Precision in the Visual Pathway through the Interplay of Excitation and Stimulus-Driven Suppression[END_REF] for some γ > 0 and even when ξ > 0.

Otherwise for other choices of nonlinearity we can get some information about the location of the solution (see Table 2).

Solving problem (16b)

Let us remind sightly rewrite problem (16b) as

u (k+1) = arg min u E α (z (k+1) , u) + 1 2γ (k) u -u (k) 2 2 , = arg min u α 2 s× d u -z (k+1) 2 2 + 1 2γ (k) u -u (k) 2 2 + u 1 + µ Hu 1 , (19)
where parameters α , µ and γ (k) are normalized by λ. The difficulty in [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF] comes from the prior term which the sum of a non differentiable separable term (the sparsity constraint, u 1 ), and non differentiable and non separable term (the smoothness constraint, µ Hu 1 ). Because of this non differentiability, one needs algorithms taken from non-smooth convex optimization.

In this paper we use a proximal algorithm described in Sect. 2.6.1.

We give here two ways to compute an approximation of (16b). The first is based on a regularization of the smoothness term and needs to introduce a small positive parameter ε but it is more computationally fast (Sect. 2.6.2). The second relies on duality and is largely inspired from algorithm used in image processing [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF][START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF] (Sect. 2.6.3).

Proximal algorithm: Forward-backward algorithm with splitting

Proximal algorithms generalize descent algorithm with projection. For example, in the particular case of an indicator function of a convex set C, the proximal operator reduces to the projector on the set C denoted P C . Generally, to solve min C g = min(g 1 + i C ) where g 1 is a differentiable function and i C the characteristic function associated to C, we use the iteration

x (n+1) = P C (x (n) -δ n ∇g 1 (x n )).
The idea of the proximal operator is to handle a more general problem by replacing the characteristic function by a general proper function g 2 and the projector P C by the proximal mapping prox γg2 with γ small enough. In [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] a collection of proximal splitting methods is presented.

Among this class of proximal algorithms, we focused on the so-called forward-backward algorithm which consists in minimizing a non differentiable convex and proper function g which decomposes as

g = g 1 + g 2 , (20)
with g 1 differentiable with a L-Lipschitz continuous gradient on its domain and g 2 a simple ‡ convex proper function. It can be decomposed in two steps: a forward (explicit) gradient step using the function g 1 and a backward (implicit) step using the function g 2 . The forward-backward ‡ We say that a function is simple if its proximal operator can be easily computed.

RR n°8837

algorithm can be interpreted as a combination of the gradient method scheme x (n+1) = x (n)δ n ∇g 1 (x (n) ) and of the proximal point iteration x (n+1) = prox δng2 (x (n) ).

Another point of view can be obtained by discretizing and searching stationary points of the dynamical system dx dt = -(∂g 1 (x) + ∇g 2 (x)). If we take an implicit scheme in ∂f 2 and an explicit one in ∇g 1 we obtain

x (n+1) = x (n) -δ n ∂g 1 (x (n+1) ) + ∇g 2 (x (n) ) ,
and from (6) this rewrites as

x (n+1) = prox δng2 (x (n) -δ n ∇g 1 (x (n) ). (21)
Algorithm ( 21) is known to converge if g 1 ∈ C 1,1 with gradient L-Lipschitz and δ n ≤ 1 L . After this short recalling on proximal algorithms, we present a variant introduced by Nesterov and revised by Beck and Teboule called Fast Iterate Shrinkage-Thresholding Algorithm (FISTA) [START_REF] Beck | Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems[END_REF]:

Algorithm 1 FISTA algorithm [2]. 1. Initialization: x 0 ∈ dom(g 1 ), y (1) = x 0 , t 1 = 1, 0 < h ≤ 1 L 2. x (n) = prox hg2 y (n) -h∇g 1 (y (n) ) 3. t n+1 = tn+ √ t 2 n +4t 2 n 2 4. y (n+1) = x (n) + tn-1 tn+1 (x (n) -x (n-1) )
Algorithm 1 will be used in the two subsequent sections. The interest of Algorithm 1 compared to the original original formulation ( 21) is convergence speed. Algorithm ( 21) is known to converge in terms of the objective function in O(1/n) [START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF]. Nesterov showed that this speed can be improved and reach O(1/n 2 ) by choosing an appropriate steps size (δ n ) n and a weighted average on the previous iterate x (n) and the candidate given by [START_REF] Gerstner | Spiking Neuron Models[END_REF]. Besides Nesterov demonstrated that it was not possible to obtain a better convergence speed by using a first order algorithm (i.e., which only use the first order derivatives of the function). More precisely, if x is the minimum of g, given an initial data x 0 ∈ dom(g 1 ) algorithms defined by Nesterov [START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF] ensure that

g(x (n) ) -g(x ) ≤ L x -x 0 2 n 2 .
We refer the reader to [START_REF] Su | A differential equation for modeling nesterov's accelerated gradient method: Theory and insights[END_REF] for a justification of the convergence of this algorithm by using dynamical system tools.

Solving problem (16b) by approximation of the smoothness constraint

As stated before, the non smooth part should be simple in the sense that its proximal operator is easily computable. On that point, since Hu 1 is not separable, it is difficult to compute its proximal operator while for the term u 1 it is easy to do. A solution is to introduce a parametrized regularization Hu 1,ε with ε > 0 small enough such that

Hu 1,ε = i,j,k ε 2 + |(Hu) i,j,k | 2
Problem ( 19) is then replaced by

min u∈X g(u) := α 2 s× d u -z 2 2 + µ Hu 1,ε + u 1 + 1 2γ (k) u -u (k) 2 2 . ( 22
)

Inria

To use Algorithm 1, let us split g as in [START_REF] Enroth-Cugell | The contrast sensitivity of retinal ganglion cells of the cat[END_REF], with

   g 1 (u) = α 2 s × d u -z (k+1) 2 2 + µ Hu 1,ε + 1 2γ (k) u -u (k) 2 2 , g 2 (u) = u 1 , (23)
where g 1 is differentiable on X and its gradient is given by

∇g 1 (u) = α s × d (s × d u -z (k+1) ) + µ H T Hu ε 2 + |Hu| 2 + 1 γ (k) (u -u (k) ), (24)
where

×
L = 144µ ε + α s 2 1 + 1 γ (k) . ( 25
)
The proximal function associated to γg 2 for some constant γ > 0 is defined by the element-wise soft threshold operator:

prox γg2 (u) = prox γ . 1 (u) = (prox γ|.| (u 1 ), ..., prox γ|.| (u N )),
where N denotes the number of elements in u and with prox γ|.| defined by (see Fig. 4)

prox γ|.| (x) =      x -γ if x ≥ γ, x + γ if x ≤ -γ, 0 otherwise. ( 26
) prox |.| (x)
x φ((M u) i ) with φ a C 1,1 convex function whose its gradient is Lispschitz, and M a matrix, we can apply Algorithm 1 by putting this term into the regular part (g 1 ).

Exact solving of problem (16b)

One problem of the previous approach is that it needs another parameter ε that cannot be fixed a priori before knowing the range of values of the unknown. To avoid this, we propose another standard approach in non-smooth optimization which consists in introducing a dual variable using the Fenchel transform [START_REF] Bolte | The Łojasiewicz Inequality for Nonsmooth Subanalytic Functions with Applications to Subgradient Dynamical Systems[END_REF].

Let us first rewrite problem [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF] as

u (k+1) = arg min u 1 2 A (k) u, u -u, b (k) + Ku µ , (27)
with K = I ×H mapping X into X ×X ν with ν = 9 for the Hessian operator in 3D (ν = 3 for the gradient in 3D),

v µ = v 1 1 + µ v 2 1 , Au = α s × d s × d u + u γ (k) and b = α s × d z (k+1) + u (k) γ (k) .
To shorten notations, we will skip the dependence on the parameter k in the sequel. We set ρ > 0 such that ρ A 2 < 1. Thanks to Lemma 2.1 we get that ( 27) rewrites as

inf u,w 1 2ρ u -w 2 2 + Cw, w -u, b + Ku µ ,
where C = (I -B) -1 B with B = ρA. We define

F (u, w) = H(u, w) + Ku µ with H(u, w) = 1 2ρ u -w 2 2 + Cw, w -u, b .
Thanks to the convexity of F (u, w), we have the relation (u, w) minimizes F (u, w) ⇔ u minimizes F (., w), w minimizes F (u, .).

Thanks to Lemma 2.1, w u = (I -B)u minimizes F (u, .). Now let us minimize F (w, .) at fixed w ∈ X . We define the two following convex sets:

H K = {K T r, r ∈ H} with H = {(p, q) ∈ X × X ν , |p i | ≤ 1, |q i | ≤ µ } . Lemma 2.2. The problem inf u 1 2ρ u -c 2 2 + Ku µ , (28)
has the unique solution u = (I -Π ρH K )(c) where Π H K is the orthogonal projection on the set

H K = {K T p, p ∈ H}. Moreover let v be such that K T v = 1 ρ Π ρH K (w), then v is solution of inf v∈H ρK T v -c 2 2 . ( 29
)
Proof. By writing the dual problem of (28) and by using relation [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF], we get:

inf u 1 2ρ u -c 2 2 + δ H K (u),
where δ H K is the support function given by [START_REF] Bolte | Clarke Subgradients of Stratifiable Functions[END_REF]. Interchanging sup and inf by using minimax theorem [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF] (VI Prop. 2.1), we get the result. For more details see [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF].

Inria

Thanks to Lemma 2.2, we deduce that u w = (I -Π ρH K )(w + ρb) minimizes F (., w). Now we can derive the following algorithm to compute the solution of (19): Proposition 2.4. Problem (27) can be solved with the following iterative algorithm:

c n = (I -B)u n + ρb, (30a)
v n = arg min v∈H c n -ρK T v 2 2 = Π ρH K (c n ), (30b)
u n+1 = c n -ρK T v n . (30c)
Remark that to compute the projection Π H ρK one needs to find the solution of ( 29) and for that we use Algorithm 1 with the expression of the projector Π H : X × X ν → H given by

(Π H ((p, q))) i = p i max(1, |p i |) , µ q i max(µ , |q i |) . ( 31
) Theorem 2.2. If ρ is such that ρ A T A 2 < 1,
then all the sequence (u n ) n converge to the solution of [START_REF] Ekeland | Analyse convexe et problèmes variationnels[END_REF].

Proof. Algorithm (30) can be rewritten as

u n+1 = (I -Π ρH K ) ((I -B)u n + ρb), (32)
which is a contraction on X (see [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF] for more details).

Conclusion: the complete algorithm

In this section we recaptulate the complete algorithm we use to solve problem [START_REF] Carandini | Do we know what the early visual system does[END_REF]. It is presented in Algo. 2. The block dedicated to the computation of u (k+1) can be changed using Sect. 2.6.2 (introducing another parameter > 0 close to 0 to regularize the term Hu 1 ) but all results presented here are obtained from Algo. 2.

Numerical results

In this section we show results using a simulated spike train, that is obtained thanks to an artificial neuron following here the LNP model (P). Section 3.1 presents two approaches from the stateof-the-art that we will use for comparison. Section 3.2 describes the experimental protocol and in particular how simulated spike train is generated. Section 3.3 gives results concerning both estimation quality and convergence.

Comparison with the state-of-the-art

In this section we remind two state-of-the-art methods used in Sect. 3.3 to make comparisons with our variational nonconvex approach (referred by VarNCvx in this section). § With g 1 (v) = ρK T v -c 2 2 and g 2 (v) = i H (v) (of proximal operator given by ( 31)).

RR n°8837

Algorithm 2 Final algorithm.

Require: Parameters of the functional ( 14): α, λ, µ, f verifying Hypotheses 1 and 2(i).

Require: Parameters of the algorithm:

N alt , N nest , N proj , (γ (k) ) 1≤k≤N alt , (β (k) ) 1≤k≤N alt , Require: Initialization: u (0) ∈ X for k = 1 : N alt do for i = 1 : N t do z (k+1) i
defined by [START_REF] Demengel | Fonctions à hessien borné[END_REF] and computed thanks to Table 1 and2

. end for Set u 0 = u (k) for n = 1 : N nest do Compute c n (30a) Compute v n (30b) with N proj iterations of Algorithm 1 § Compute u n (30c) end for Set u (k+1) = u Nnest end for

Variational approach with convex nonlinear function f

As explained in the introduction which has motivated all this work, in state-of-the-art variational approaches, convex nonlinear functions f are used [START_REF] Paninski | Maximum likelihood estimation of cascade point-process neural encoding models[END_REF][START_REF] Mcfarland | Inferring Nonlinear Neuronal Computation Based on Physiologically Plausible Inputs[END_REF][START_REF] Pillow | Likelihood-based approaches to modeling the neural code. Bayesian brain: Probabilistic approaches to neural coding[END_REF]. To show what is the interest of using a non convex function f , the same Algorithm (16a)-(16b) will be applied but with a convex function instead. This method will be referred to as VarCvx and based on the following function f c (see Fig. 5(e)):

f c (x) = 0 if x < -1 2 , 1 2 + 2x + 2x 2 otherwise. (33)

Spike Triggered average (STA)

The STA approach is a very classical method used in the neuroscience community to obtain a discrete approximation of the receptive field of sensory neurons [START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF].

In the case of a visual neuron stimulated by a visual stimulus s, STA estimator is defined by

STA(x, τ ) = 1 n(T ) n(T ) i=1 s(x, t i -τ ),
where {t i } 1≤i≤n(T ) is the sequence of n(T ) spikes times generated by the neuron.

A common visual stimulus used in experiments is a sequence of white noise images (see Fig. 5(a) and also Fig. 1 for illustrations), i.e., images whose the power spectrum is constant. That particular case is interesting since it can be shown that the STA estimator at time n(T ) fixed converges in law to N (u, σ 2 n(T ) ) when n(T ) is large and where u is the receptive field to recover. Even if the process describing the stimulus does not generate exactly a white noise but only an independent discrete uniform binary random signal, by an heuristic way, thanks to the central limit and Bussgang's theorems, convergence remains true when n(T ) → +∞. Remark 3.1. An important parameter of the STA is the spatial resolution of the estimation, which corresponds to the size of the blocks of the white noise stimulus images. In general, it is experimentally fixed to reach a compromise to have a neuron sufficiently responding: indeed, for Inria STA to be precise, one needs a sufficient number of spikes, which can be obtained both by increasing the experiment duration (still limited by physiological constrains) or choosing an "optimal" block size: If too small, neuronal responses might be too weak since receptive field would have to integrate too small details, thus leading to receptive field with low Signal-to-Noise-Ratio; on the contrary, if too large, smallest receptive field will be lost, or not described with enough details, because of the coarse approximation.

Simulated spike train: the direct problem

In this section we describe how simulated spike train is generated, following the direct problem described in Fig. 1 and5 and assuming that we consider classical ON type retinal ganglion cell which can be approximated by a LNP model. ¶ Stimulus We choose as visual stimulus a sequence of N t binary white noise images of size 20 × 20 pixels with block size 4 × 4 pixels, thus allowing direct comparisons with STA in the sequel. One sample image is shown in Fig. 5(a). Each images is fed as input to the neuron during a time period of length ∆t. Following notations of Sect. 2, up to changing the maximum of the sigmoid c by c∆t, we choose ∆t = 1.

Receptive field of the artificial neuron

The linear part of the response of classical retinal ganglion cells is generally modeled by a nonseparable spatio-temporal filtering (see, e.g., [START_REF] Wohrer | Virtual Retina: A biological retina model and simulator, with contrast gain control[END_REF]). A common approximation is to consider that it is separable in space and time. We make this assumption here. Note that this is only to define our ground truth more simply, but of course this is not an hypothesis needed for our approach to work.

We assume that the receptive field u of the ON retinal ganglion cell is separable in space and time, so that u(x, t) = v(x)w(t).

The spatial part can be described by a difference of Gaussian functions (DoG) (Fig. 5(b)-(c)):

v(x) = ω C G σ C (x) -ω S G σ S (x),
where

G σ (x) = e -|x| 2 2σ 2 .
This center-surround behavior which has been found in the 60's [START_REF] Rodieck | Quantitative analysis of cat retinal ganglion cell response to visual stimuli[END_REF][START_REF] Enroth-Cugell | The contrast sensitivity of retinal ganglion cells of the cat[END_REF] corresponds to a measure of the local contrast. This DoG is driven by four parameters which can be all interpreted functionally. σ C is a measure of the blur applied to the image hitting the retina originating in particular from the sampling frequency of photoreceptors. In some sense it defines the precision of the retina. Ratio σ S /σ C defines the relative surround. ω C is a linear gain which gives the orders of magnitude for retina amplification. Ratio ω S /ω C defines the relative surround weight.

In our simulations, we chose σ C = 2.2, σ S = 3, ω C = 100 and ω S /ω C = 0.9.

The temporal part can be described by a difference of Exponential functions (DoE) (Fig. 5(d)):

w(t) = ω A E(t, n A , τ A ) -ω B E(t, n B , τ B ), with τ A < τ B , ¶
Remark that primate retina has around 20 types of ganglion cells which interrogates any point in visual space about a number of distinct qualities [START_REF] Wohrer | Model and large-scale simulator of a biological retina with contrast gain control[END_REF][START_REF] Gollisch | Eye smarter than scientists believed: neural computations in circuits of the retina[END_REF][START_REF] Masland | Cell populations of the retina: The proctor lecture[END_REF]. This is achieved via parallel complex retinal circuitry which cannot be all approximated accurately by LNP models.

We choose as convention that s takes values in {-1, 1}.

where

E(t, n, τ ) = (nt) n exp(-nt τ ) (n -1)!τ n+1 .
This band-pass temporal filter has a strong functional similarity with the DoG model. The first part E(t, n A , τ A ) corresponds to the low-pass properties while the second part E(t, n B , τ B ) represents the delayed inhibition that makes the response transient. * * Ratio ω B /ω A defines the strength of the transient. In our simulations we chose

ω A = ω B = 1, τ A = 5, τ B = 7, n A = 5, n B = 7
Receptive field is then discretized as a spatio-temporal volume of size 20 × 20 × 30 pixels. It defines our ground truth (GT). Some temporal slices are shown in Fig. 6 (first row). So, given u, one can already compute the linear response s × d u.

Nonlinearity of the artificial neuron Given s × d u, one has to set the nonlinearity f to compute the rate of the Poisson process generating spikes (P). To do so, we choose to take f (x) = cf 0 (ax + b) where f 0 is the piecewise cubic function [START_REF] Su | A differential equation for modeling nesterov's accelerated gradient method: Theory and insights[END_REF], so that the distribution of s × d u values mainly falls outside saturation regime, so that neuron response present sufficient variability to recover receptive field (see Fig.

Reconstruction

This section presents results obtained with our approach, using Algorithm 2 described in Sect. 2.7. Parameters of the approach were chosen as α = 1000, γ (k) = β (k) = 10, N alt = 300, N nest = 40, * * In response to an input step function, DoE temporal filtering part makes induces a peal of activity, and then a decrease of activity as the inhibition (of time scale τ B ) builds up and compete with the direct signal (of time scale τ A ).

Inria

N proj = 20 and u 0 = 0.

Figure 6, shows comparisons of receptive field estimation with different approaches. Second row shows STA estimation. Result is noisy and resolution is constrained by the block size of the stimulus (by definition). Remark that the block size of the stimulus should depend on the receptive field to reconstruct (in general, it cannot be of size 1 pixel) because of the operator × d which sums spatially on the stimuli. Hence if the support of the receptive field is composed of several pixels, the linear response s × d u leads to a gray value constant in time, so that receptive field cannot be recovered (see also Remark 3.1). Third row shows a linear interpolation of receptive field to have a 20 × 20 pixels spatial resolution. This is given to compare with results from variational approaches which have that resolution. Fourth and fifth rows show results using the variational approach, with the convex nonlinearity f c (VarCvx) and the nonconvex nonlinearity f (VarNCvx) respectively. Qualitatively, variational approach enables to reconstruct very accurately the receptive field compared to STA. Note that significant qualitative differences between VarCvx and VarNCvx appear. In particular, VarCvx does not enable to reconstruct the surrond of the receptive field while VarNCvx does (see Fig. 8 for a more quantitative comparison).

Figure 7 shows the influence of weights λ and µ defining the prior in energy [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF]. Qualitatively, as expected, the more the sparsity weight λ is large the more the receptive field is localized, and the more the smoothness constraint is large, the more the receptive field is smooth. In both cases, we find one optimal value of parameters to reach the lowest covariance error.

Figure 8 compares the temporal profiles of receptive field obtained with VarCvx and VarNCvx. Results show that VarCvx case is always biased while VarNCvx converges to the ground truth when the duration of the experiment increases (equivalently when the number of spikes increases).

Figure 9(a)-(b) compare the rate of convergence w.r.t the number of spikes observed. Our approach converges faster in the sense of both the covariance error E cov (u, u GT ) and the l 2 -norm E l2 (u, u GT ). The three methods converge in the sense of the covariance error but we see that in the sense of the l 2 error the variational approach with the convex approach does not converge. Another observation is that error decays linearly (in log scale) for both the STA approach and the nonconvex variational approach with a slope of the convergence equal for both methods which can be explained by the central limit Theorem.

Figure 9(c) illustrates the speed of convergence of the algorithm: we compare the error decaying for our approach VarNCvx and for VarCvx, using a convex approximation f c (33) as usually done in the literature.

Finally, Fig. 10 shows the execution time in seconds of Algorithm 2 implemented in Matlab and run on a GPU Quadro K4000.

Conclusion

Following the classical model of visual sensory neuron that is the linear-nonlinear Poissonian (LNP) model, we focused on the case when the nonlinearity is a sigmoid which is more realistic from a physiological point of view. This work presents for the first time, up to our knowledge, an approach to solve the nonconvex variational formulation for receptive fields estimation. Compared to the state-of-the-art, results on synthetic data are very promising. Our approach allows to estimate receptive field with great accuracy and it converges faster w.r.t. stimulus duration, which makes it a good candidate for experimentalists who need to estimate receptive fields in a limited amount of time † † . GT

µ E cov (u, u GT ) E cov (u, u GT ) (a) (b)
new (casbis, with stronger surround)

Figure 7: Quantitative influence of prior term parameters on receptive field estimation, using VarNCvx. We show covariance error (E cov (u, u GT ) = 1cov(u, u GT )/σ u σ u GT ) as a function of (a) the sparsity constraint parameter λ (for µ = 1 fixed) and (b) the regularity constraint parameter µ (for λ = 1 fixed).

E cov (u, u GT ) n(T ) n(T ) (log scale) E l2 (u, u GT ) (log scale) E l2 (u, u GT )
E cov (u, u GT ) n(T ) n(T ) (log scale) E l2 (u, u GT ) (log scale) E l2 (u, u GT )

RR n°8837

Future work will focus on validating our approach on real cell recordings. Our plan is to study recordings coming from a new generation of large-scale, high density MEA consisting of 4096 electrodes (APS CMOS-MEA, [START_REF] Berdondini | Active pixel sensor array for high spatio-temporal resolution electrophysiological recordings from single cell to large scale neuronal networks[END_REF][START_REF] Maccione | Following the ontogeny of retinal waves: pan-retinal recordings of population dynamics in the neonatal mouse[END_REF]. The main interest of this technology will be to provide recordings of thousands of cells, allowing to benchmark our approach w.r.t. the state of the art based on a large population of cells.

This work also opens new avenues to further improve it. Our approach works for any kind of stimulus, contrarily to classical approaches such as STA which assume that stimulus should be sequences of white noise images. So this allows to explore and optimize stimulus itself to improve reconstructions and convergence (see, e.g., [START_REF] Pamplona | A super-resolution approach for receptive fields estimation of neuronal ensembles[END_REF][START_REF] Hilgen | Pan-retinal characterization of light responses from ganglion cells in the developing mouse retina[END_REF] for an example of alternative stimulus used to increase resolution). Here, we focused on LNP model, however, it is well known that some sub-population of cells cannot be modelled with such model (e.g., ON-OFF cells which respond to both onset and offset of a stimulus). So it will be interesting to investigate how our formulation can be extended using models generalizing LNP, e.g., by adding a non-linear feedback on spikes [START_REF] Pillow | Spatio-temporal correlations and visual signalling in a complete neuronal population[END_REF], adding a suppresive LNP cascade [START_REF] Butts | Temporal Precision in the Visual Pathway through the Interplay of Excitation and Stimulus-Driven Suppression[END_REF][START_REF] Mcfarland | Inferring nonlinear neuronal computation based on physiologically plausible inputs[END_REF]. Finally, when applying our approach to large populations of cells, one may need to have faster execution time which should be possible thanks to stochastic algorithm with asychronous parallelization [START_REF] Peng | ARock: an Algorithmic Framework for Async-Parallel Coordinate Updates[END_REF].

Proof. Let dt small enough, we have by denoting g k (t) = P (n(t)n(t 0 ) = k):

g n (t) = n k=0 g n-k (t)P (n(t + dt) -n(t) = k) = g n (t)P (n(t + dt) -n(t) = 0) + g n-1 (t)P (n(t + dt) -n(t) = 1) + R n (t, dt) = g n (t)(1 -λ(t)dt + o(dt)) + g n-1 (t) (λ(t)dt + o(dt)) + R n (t, dt), (35)
with

R n (t, dt) = n k=2 g n-k (t)P (n(t + dt) -n(t) = k) = o(dt).
By passing to the limit in [START_REF] Pillow | Likelihood-based approaches to modeling the neural code. Bayesian brain: Probabilistic approaches to neural coding[END_REF] when dt → 0 we get

dg n dt + λg n = λg n-1 . (36)
By using that g 0 (0) = 1, the solution of this differential equation is

g n (t) = Λ(t 0 , t) n n! e -Λ(t0,t) ,
with Λ(t 0 , t) = t t0 λ(s)ds.

Proof of formula (3). Let us first define a discretized approximation of u(x, t). To do so, we split Ω × [0, T [ in N small voxels and we denote u N i the constant approximation of u(x, t) in the i-th voxel. Similarly, we split [0, T [ in N t small intervals of length ∆t = T Nt and we denote by ξ Nt (t) the function equal to ξ Nt i for t in [(i -1)∆t, i∆t[, where ξ Nt i is the number of spikes in that interval. We have that

ξ Nt (t) -→ Nt→∞ n(T ) i=1 1 t=ti , u N (x, t) -→ N →∞ u(x, t),
in the simple sense. Now we denote by ξ Nt i the value of ξ Nt (t) when t belongs to the i-th bin and by u N i the value of u N (x, t) on the i-th voxel. We denote by u N , s N λ Nt and ξ Nt the vectors of component respectively

u N i , s N i , λ Nt i = f ([s N × d u N ] i
) and ξ Nt i . From Lemma A.1 and independence of ξ Nt i , we have

P (ξ Nt |u N ) = Nt i=1 P (ξ Nt i |u N ) = e -∆t k λ N t k Nt i=1 (λ Nt i ∆t) ξ N t i ξ Nt i !
.

Let be N t tending to infinity, we get that

ρ({t i } 1≤i≤n(T ) |u N ) = e -T 0 λ N n(T ) i=1 λ N (t i ),
where λ N (t) = f ([s × u N ](t)). Let be N → ∞ we get the well known formula [START_REF] Paninski | Maximum likelihood estimation of cascade point-process neural encoding models[END_REF]. Passing to the -log we get

-log(ρ({t i } 1≤i≤n(T ) |u)) = T 0 λ(s)ds - n(T ) i=1 λ(t i ). (37)
B Study of the operator ×

In this appendix we study properties of the operator × in the continuous case. Note that same estimations can be obtained in discrete case. To simplify notations, up to extend the temporal functions by 0 outside their support, we assume that the temporal support is R and we denote by Ω ⊂ R n the spatial domain. Notations which are used here are specific to this appendix (f denotes whatever function and is not linked the nonlinearity used in our model (P)). The following lemma is an adaptation of the Young inequality proof (see [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF], Th. 4.15) to the infinite dimensional case (here the state space is L q (Ω) for some q ≥ 1).

Lemma B.1. Let be s ∈ L 1 (R, L q (Ω)) and u ∈ L 1 (R, L q (Ω)) with q such that 1 q + 1 q = 1, then

s × u L 1 (R) ≤ s L 1 (R,L q (Ω)) u L 1 (R,L q (Ω)) .
Proof. We set F (t, τ ) = Ω s(x, tτ )u(x, τ )dx, by Holder inequality we have

|F (t, τ )| ≤ s(t -τ ) L q (Ω) u(τ ) L q (Ω) .
By summing over t we get

R |F (t, τ )|dt ≤ s L 1 (R,L q (Ω)) u(τ ) L q (Ω) ,
and by summing over τ we obtain the result.

Corollary B.1. Let s ∈ L 1 (R, L q (Ω)) and u ∈ L p (R, L q (Ω)) with q such that 1 q + 1 q = 1 and

1 ≤ p ≤ ∞, then s × u L p (R) ≤ s L 1 (R,L q (Ω)) u L p (R,L q (Ω)) .
Proof. Let us assume that 1 < p < ∞ and let p be such that 1 p + 1 p = 1. Keeping the same notations as in the proof of Lemma B.1, one can rewrite the upper bound of |F (t, τ )| as

|F (t, τ )| ≤ s(t -τ ) 1 p L q (Ω) ∈L p w.r.t. τ s(t -τ ) 1 p L q (Ω) u(τ ) L q (Ω) ∈L p w.r.t. τ .
By summing over τ and thanks to Holder inequality we have

R |F (t, τ )dτ | ≤ s 1 p L 1 (R,L q (Ω)) s L q (Ω) u p L q (Ω) (t) 1 p ,
where denotes the usual convolution operator. Taking the power p we get

|(s × u)(t)| p ≤ s p p L 1 (R,L q (Ω)) s L q (Ω) u p L q (Ω) (t).
By using Lemma B.1 we have

s L q (Ω) u p L q (Ω) L 1 (R) ≤ s L 1 (R,L q (Ω)) |u| p L 1 (R,L q (Ω)) , = s L 1 (R,L q (Ω)) u p L p (R,L q (Ω)) .
Summing over t and taking the p-th root lead to the result.

Inria

Lemma B.2. Let s ∈ L p (R, L q (Ω)) and z ∈ L p (R) with 1 p + 1 p = 1, then the adjoint operator of × is given by

(s × z)(x, t) = R s(x, τ -t)z(τ ),
and s × z ∈ L p (R, L q (Ω)) with the estimation

s × z L p (R,L q (Ω)) ≤ s L 1 (R,L q (Ω)) z L p (R) .
Proof. Computing the scalar product of s × u and z gives the expression of × . Then the estimation of the norm is similar to the proof of Cor. B.1.

Remark B.1. When q = p = 2 and the sequence s ∈ L 1 (R, L 2 (Ω)) is finite, from Lemma B.2 and Corollary B.1, the Lipschitz constant of s × . and s × . in norm L 2 (R) and

L 2 (R, L 2 (Ω)) is bounded by s L 1 (R,L 2 (Ω)) . Hence the Lipschitz constant of the operator g → s × s × g from L 2 (Ω) into itself is bounded by s 2 L 1 (R,L 2 (Ω)) .

C Around the Kurdyla-Lojasiewicz property

To prove convergence, we apply Th. 9 from [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF]. The difficult point is to check that E α verifies the Kurdyka-Lojasiewicz (KL) property which is a notion introduced in 1963 by Lojasiewicz for analytical functions and extended by Kurdyka in 1998 to functions definable in an o-minimal structure. To be self-content, let us first recall some definitions and notions associated to this property.

Definition C.1 (KL property). We say that a proper function φ has the Kurdyka-Lojasiewicz (KL) property at x ∈ dom(∂φ) if there exist a neighborhood V of x, ν ∈]0, ∞] and a continuous concave function ψ : [0, ν[→ R + with ψ(0) = 0 and such that:

(i) ψ is continuously differentiable on (0, ν) with ψ > 0.

(ii) For all x ∈ V with φ(x) ≤ φ(x) ≤ φ(x) + ν, one has

ψ (φ(x) -φ(x))dist(0, ∂φ(x)) ≥ 1. (KL)
A proper closed function f satisfying the KL property at all points in dom(∂φ) is called a KL function.

Definition C.2 (KL exponent). For a proper function φ satisfying the KL property at x ∈ dom(φ), if the corresponding function ψ can be chosen as ψ(s) = cs 1-α for some c > 0 and α ∈ [0, 1[, i.e., there exists c, ε > 0 and ν ∈]0, ∞] such that dist(0, ∂φ(x)) ≥ c(φ(x)φ(x)) α , whenever xx 2 ≤ ε and φ(x) ≤ φ(x) ≤ φ(x) + ν, then we say that φ has the KL property at x with an exponent α. If φ is a KL function and has the same exponent α at x ∈ dom(∂φ), then we say that φ is a KL function with an exponent of α.

Remark C.1. A proper lower semicontinuous function φ : R n → R∪{+∞} has the KL property at every non critical point (see Lemma 2 in [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF]) for whatever exponent α ∈ [0, 1].

The exponent of the function involved in the KL inequality gives information on the rate of convergence toward critical points for descent algorithms used in optimisation [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF] and toward stationary points for subgradient trajectories [START_REF] Bolte | The Łojasiewicz Inequality for Nonsmooth Subanalytic Functions with Applications to Subgradient Dynamical Systems[END_REF] from a dynamical systems point of view.

In order to determine which functions satisfy the KL property we introduce the notion of semialgebraic sets. A subset of R n is called semialgebraic if it can be written as a finite union of sets of the form {x ∈ R n , p i (x) = 0, q i (x) < 0, i = 1, ..., p}, } where p i , q i are real polynomials. A proper function φ : R n → R ∪ {+∞} is semialgebraic if its graph is a semialgebraic subset of R n+1 . It is known that semialgebraic functions verify the KL property.

We recall the definition of o-minimal structure introduced in [START_REF] Dries | Tame topology and o-minimal structures[END_REF] which can be seen as an axiomatization of properties of semialgebraic sets.

Definition C.3. Let O = {O n } n∈N be such that each O n is a collection of subsets of R n . The family O is an o-minimal structure over R if it satisfies the following axioms:

(i) Each O n is a boolean algebra. Namely ∅ ∈ O n and for each A, B in O n , A ∪ B, A ∩ B and R n \A belong to O n . (ii) For all A in O n , A × R and R × A belong to O n+1 . (iii) For all A in O n+1 , Π(A) := {(x 1 , ..., x n , ) ∈ R n , (x 1 , ..., x n , x n+1 ) ∈ A} belongs to O n . (iv) For all i = j in {1, ..., n}, {(x 1 , ..., x n ) ∈ R n , x i = x j } ∈ O n . (v) The set {(x 1 , x 2 ) ∈ R 2 , x 1 < x 2 } belongs to O n .
(vi) The elements of O 1 are exactly finite unions of intervals.

Given an o-minimal structure O, a set A is said to be definable (in O), if A belongs to O. A function φ : R n -→R ∪ {+∞} is said to be definable if its graph is a definable subset of R n × R.

o-minimal structures have interesting properties:

• Finite sums of definable functions are definable.

• Indicator functions of definable sets are definable.

• Compositions of definable functions are definable.

• Generalized inverses of definable functions are definable.

• Functions of the type

R n x → f (x) = sup y∈C g(x, y) or R n x → f (x) = inf y∈C g(x, y)
where g and C are definable, are definable.

The class of semialgebraic sets is an o-minimal structure [START_REF] Benedetti | Real algebraic and semi-algebraic sets[END_REF] that we denote semi-alg.

A subset X ⊂ R n is semianalytic if at each point a ∈ X there exists a neighborhood U such that X ∩ U can be written as finite union and intersection of analytic equalities and inequalities. A subset X ⊂ R n is subanalytic if each point of R n admits a neighborhood U such that X ∩ U can be written as the projection of a bounded semianalytic subset of R n × R m for some m ≥ 1.

Notice that image and preimage of a subanalytic set are not in general a subanalytic set. Let τ n defined by

τ n (x 1 , ..., x n ) = x 1 1 + x 2 1 , ..., x n 1 + x 2 n Inria A subset X ⊂ R n is globally subanalytic if its image under τ n is a subanalytic subset of R n .
Globally subanalytic sets are subanalytic and reciprocally any bounded subanalytic sets are globally subanalytic. The class of globally subanalytic sets is an o-minimal structure that we denote global-subanal. There exists an o-minimal structure denoted log-exp [START_REF] Wilkie | Model completeness results for expansions of the ordered field of real numbers by restricted Pfaffian functions and the exponential function[END_REF][START_REF] Dries | Tame topology and o-minimal structures[END_REF] which contains global-subanal and the graph of exp : R → R. This huge o-minimal structure contains all the aforementioned structures.

The following theorem enables to characterize KL functions by using the notion of o-minimal structure.

Theorem C.1 (from [START_REF] Bolte | Clarke Subgradients of Stratifiable Functions[END_REF]). Any proper lower semicontinuous function φ : R n → R ∪ {+∞} which is definable in an o-minimal structure O has the KL property at each point of dom(∂φ). Moreover the function ψ of the (KL) inequality is definable in O.

Proof of Theorem 2.1

In order to apply Th. 9 of [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF], the only point to check is if E α verifies the KL property (see Def. C.1). Thanks to Th. C.1, a sufficient condition is that E α be definable in an o-minimal structure. Hence we need to find the smallest structure in which our functional is defined. Our energy is composed of three kinds of terms:

ϕ 1 (z) = ψ ξ (z), ϕ 2 (z, u) = z -s × d u 2 2 , ϕ 3 (u) = λ u 1 + µ Hu 1 .
For each term, one can show that it is definable:

• With Hypothesis 2, in the case ξ = 0 the first term ϕ 1 (z) is definable in the log-exp structure. For ξ > 0 by Hypothesis 2 and stability by composition with log and linear transform, ϕ 1 (z) is still be definable in this structure.

• The second term is semialgebraic so it is definable in the log-exp structure.

• The third term is the sum of terms of the form |(Hu

) i | R ν where the norm |v| R ν = ν i=1 v 2 i
is semialgebraic (its graph can be written as

{(v, z), ν i=1 v 2 i = z 2 , z ≥ 0}.
Hence thanks to stability by sum and composition with linear operator, we deduce that ϕ 3 (u) is definable in the log-exp structure.

With the stability by the × operator (axiom (ii) of Def. C.3), we deduce that (z, u) → ϕ 1 (z) and (z, u) → ϕ 3 (u) are also in the log-exp structure. In conclusion, the energy

E α (z, u) = i ψ ξi (z) + α 2 ϕ 2 (z, u) + ϕ 3 (u)
is definable in the log-exp structure. Theorem C.1 gives that E α verifies the KL property with a function ψ in (KL) which is definable in the log-exp structure. By applying Theorem 9 of [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF] we get the result.

Remark C.2. Let us notice that to show that E α verifies KL we have shown that E α is definable in the log-exp structure. Hence we have no more information on the function ψ in (KL) than it is definable in this later structure. In general, f is piecewise defined by analytic expression, but if we can show that it admits locally an analytic extension then E α is subanalytic and verifies (KL) locally with ψ a power function [START_REF] Bolte | The Łojasiewicz Inequality for Nonsmooth Subanalytic Functions with Applications to Subgradient Dynamical Systems[END_REF]. This fact is important because according to the power of the function ψ, the rate of convergence is different (see [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF], Th. 11).

D Solving problem (16a) for a piecewise linear sigmoid

This appendix gives the expression of the solution of (16a) in the case of a piecewise linear continuous sigmoid-like non linearity (see Figure 3(b)). We will see that in this case (16a) cannot be reduced to an equality. The normalized piecewise linear sigmoid function that we study is

f (x) =              0, for x < - 1 2 , 1 2 + x, for |x| ≤ 1 2 , 1, for x > 1 2 (38)
Remark D.1. The non linearity f (38) does not verify Hypothesis 2(ii): at point x = 1 2 , f cannot be locally split into the difference of a convex function and a quadratic one. Hence the proximal operator can be multivalued (which will be the case here, see below). However, convergence of ( 16) is still guaranteed thanks to Th. 2.1.

D.1 Case ξ = 0

Let γ > 0, we are interested in computing the proximal operator associated to γf . Since f (y) is constant for |y| ≥ 1 2 , we split R into the union of two disjoint sets:

I 1 = -1 2 , 1 2
and I 2 = -∞, + 1 2 ∪ 1 2 , +∞ , and we compute:

g 1 (x) = min
y 1 (x) =      -1 2 if x ≤ -1+2γ 2 , x -γ if -1+2γ 2 ≤ x ≤ 1+2γ 2 , 1 2
if x ≥ 1+2γ 2 , and

y 2 (x) =          x if x ≤ -1 2 , -1 2 if -1 2 ≤ x ≤ γ, 1 2 if γ ≤ x ≤ 1 2 , x if x ≥ 1 2 .
Computing g(x) = min(g 1 (x), g 2 (x)) we get the expression of the proximal operator of γf (see Figure 11)

prox γf (x) =          x if x ≤ -1 2 , -1 2 if -1 2 ≤ x ≤ γ, x -γ if -1+2γ 2 + γ ≤ x ≤ 1+γ 2 , x if x ≥ 1 2 . ( 39
)
Let us underline that the proximal function of γf is multivalued at x = 1+γ 2 as we can see on Fig. 11.

D.2 Case ξ > 0

We recall the expression of the function ψ ξ [START_REF] Butts | Temporal Precision in the Visual Pathway through the Interplay of Excitation and Stimulus-Driven Suppression[END_REF] in this case, extended on {f ≤ 0} by the value +∞ (Fig. 12):

ψ ξ (x) = 1 2 + x -ξlog 1 2 + x if |x| < 1 2 , 1 if x ≥ 1 2 .
- Remark D.2. ψ ξ is convex as soon as ξ ≥ 1. This fact comes from the non increasing of ψ ξ for ξ ≥ 1.

1 2 1 2 f -1 (ξ) 1 -1 2 1 2 1 ξ < 1 ξ ≥ 1
Let γ > 0, in order to compute the proximal operator of γψ ξ , we compute where ϕ x,γ (y) = f (y)ξ i log(f (y)) + 1 2γ (xy) 2 . It is easily seen that

h 1 (x) = min
z 2 (x) = 1 2 if x ≤ 1 2 , x if x ≥ 1
2 . We set φ x,γ (y) = 1 2 + yξlog( 1 2 + y) + 1 2γ (xy) 2 for y > -1 2 then

φ x,γ (y) = 0 ⇔ y 2 + γ + 1 2 -x y + γ 2 - x 2 -γξ = 0 , ξ > 0.
We search roots of φ x,γ (y) = 0 belonging to [-1 2 , 1 2 ]. We compute ∆ = x -(γ -1 2 ) 2 + 4γξ > 0, there are two distinct roots

y ± (x) = -(γ + 1 2 -x) ± √ ∆ 2 . (40)
We check that y + ≥ -1 2 and y -≤ -1 2 for all x ∈ R. Besides

y + ≤ 1 2 ⇔ x ≤ x 1 := γ(1 -ξ) + 1 2 .
By setting x 1 = γ(1ξ) + 1 2 , we obtain the following tables of variations: y -1 We can show that x → y + (x) is not decreasing and we have y + (x 1 ) = 1 2 and lim x→-∞ y + (x) = -1 2 . We get that:

z 1 (x) = y + (x) if x ≤ x 1 , 1 2 if x ≥ x 1 ,
and by computing h = min(h 1 , h 2 ) we get the expression of prox γψ ξ (Fig. 13):

• Case ξ ≤ 1

prox γψ ξ (x) = y + (x) if x ≤ x 1 , x if x ≥ x 1 .
• Case ξ > 1

prox γψ ξ (x) =      y + if x ≤ x 1 , 1 2 if x 1 ≤ x ≤ 1 2 , x if x ≥ 1 2 .
We conclude this study by the observation that for ξ < 1 the function ψ ξ (z) is not semiconvex (even not lower C 2 ), that leads to the fact that its proximal operator is multivalued at x = x 1 and is given at this point by the set { 1 2 , x 1 }. E Solving problem (16a) for a piecewise cubic sigmoid

Inria -∞ x 1 1 2 -1 2 1 2 y 1 (x) x -∞ x 1 = 1 2 -1 2 1 2 y 1 (x) x -∞ x 1 1 2 -1 2 1 2 y 1 (x) 1 2 x ξ ≤ 1 ξ = 1 ξ > 1
The cubic approximation that we choose is:

f (x) =      0 if x ≤ -1 2 , 1 2 + 3 2 x -2x 3 if |x| ≤ 1 2 , 1 if x ≥ 1 2 . (41)
Remark E.1. The sigmoid f given by (41) verifies Hypothesis 1 and 2:

(i) f is semialgebraic (see Appendix C) since its graph can be written as a finite union of polynomial inequalities.

(ii) The function y → -log(f (y)) is convex on ] -1 2 , +∞[. (iii) The function f (y) + 3y 2 is convex.

We denote by ϕ x,γ the function to minimize at fixed x: ϕ x,γ (y) = f (y) + 1 2γ (xy) 2 , and we set g(x) = min y ϕ x,γ (y).

In practice γ will be close to 0 since α will be chosen large (see [START_REF] Demengel | Fonctions à hessien borné[END_REF]). Hence, we assume that γ ≤ 1 6 so that ϕ x,γ be convex for all x. Following the same reasoning as for the computation of the proximal operator in the case of a piecewise linear sigmoid, we compute and g(x) = min(g 1 (x), g 2 (x)) and the proximal operator is defined by prox γf (x) = y 1 (x) if g(x) = g 1 (x), y 2 (x) otherwise.

g 1 (x) = min

RR n°8837

Figure 1 :

 1 Figure 1: Illustration of the LNP model (P) (gray box). Two kinds of nonlinear functions f are illustrated: ramp-or sigmoid-like nonlinearities.Here the stimulus is a sequence of white noise images, which is the classical stimulus used for the spike-triggered averaged approach (STA, see[START_REF] Chichilnisky | A simple white noise analysis of neuronal light responses[END_REF]).

Figure 2 :

 2 Figure 2: Graphs of the data fidelity term ψ ξ (.) for ξ ∈ R + , depending on values of ξ.

Figure 3 :

 3 Figure 3: Sigmoid-like functions. (a) Exponential based function

Figure 4 : 2 . 5 .

 425 Figure 4: Graph of function prox γ|.| (x) defined by (26) which is a soft thresholding

Lemma 2 . 1 (u -w 2 2 +

 212 Bect et al. [3]). If B is a symmetric positive matrix such that B 2 < 1 then we have Bu, u = inf w Cw, w ,with C = B(I -B) -1 .Moreover, the maximum is reached at w = (I + C) -1 u = (I -B)u.

  5(e) for illustration). This is done empirically. Parameters are a = 0.167, b = 0.1, c = 0.8. Simulated spike train Since we use as temporal discretization step the inter-frame ∆t, we simulate ξ = (ξ i ) 1≤i≤Nt , the vector associated to the number of spikes occurring during each bin [(i -1)∆t, i∆t[. From Lemma A.1, ξ follows a Poisson distribution of parameter f (s × d u).

Figure 1

 1 was giving an example of simulated spike train with bars width proportional to the number of spikes by bin. Here the number of spikes is around 500.

Figure 5 :

 5 Figure 5: How simulated spike trains are generated (see Sect. 3.2)? (a) One image of the stimulus consisting in a sequence of images with white noise. This kind of stimulus is commonly used for STA estimator, and it is also used here in our approach for comparison. (b) and (c) Spatial profile of the receptive field used to generate the simulated spike train. (d) Temporal profile of the receptive field. (e) Nonlinearities super-imposed on the histogram of s × d u for two cases: there is the non convex one which is used to generate the spike train and for our approach, and the convex one used for comparion (Sect. 3.1.1).

Figure 6 :

 6 Figure 6: Qualitative comparison of receptive field reconstructions with respect to ground truth (first row, GT). Columns 1-4 show a selection of temporal slices and column 5 shows a horizontal cut passing through the center of receptive field(at temporal depth represented in column 4). Rows 2 and 3 show STA results. Rows 4 and 5 show results obtained using the variational approach, with the convex nonlinearity f c (VarCvx) and the nonconvex nonlinearity f (VarNCvx) respectively (for λ = 9 and µ = 11).

Figure 8 :

 8 Figure 8: Comparison of the temporal profile of the receptive field estimates as function of the duration of the experiment N t . Results show that VarNCvx gives a non-biaised estimated better than VarCvx.0 5 10 15 20 25 30 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12

Figure 9 :Figure 10 :

 910 Figure 9: Study of the convergence of the VarNCvx approach. (a) Covariance error (E cov (u, u GT )) as a function of the number of spikes for the convex and non convex variational approaches and the STA. (b) l 2 -norm error (E l2 (u, u GT ) = uu GT 2 ) as a function of the number of spikes in log-log scale. (c) l 2 -norm error as a function of the number of iterations for N t = 1000.

|y|≤ 1 2 ϕ 2 ϕ 2 ϕ 2 ϕ

 2222 x,γ (y) with y 1 (x) = arg min|y|≤ 1 x,γ (y), g 2 (x) = min |y|≥ 1 x,γ (y) with y 2 (x) = arg min |y|> 1 x,γ (y),where ϕ x,γ (y) = f (y) + 1 2γ (xy)2 . Standard computations lead to:

Figure 11 :

 11 Figure 11: Plots of functions used in the piecewise linear sigmoid case. (a) prox γf (x). (b)ϕ γ (x) = f (x) + 1 2γ (x -1+γ 2 ) 2 . (c) f γ (x) = inf y f (y) + 1 2γ (xy) 2 the Moreau envelop of f compared to f .

Figure 12 :

 12 Figure 12: Shape of the data fidelity term ψ ξ in the case of nonlinear sigmoid-like piecewise linear function f (38).

  (y+) case (x ≥ x 1 ) case (x ≤ x 1 )

Figure 13 :

 13 Figure 13: Plots of prox γψ ξ (x) for different values of ξ > 0.

|y|≤ 1 2 ϕ 2 ϕ 2 ϕ 2 ϕ

 2222 x,γ (y), y 1 (x) = arg min|y|≤ 1 x,γ (y), g 2 (x) = min |y|≥ 1 x,γ (y), y 2 (x) = arg min |y|≥ 1 x,γ (y),

Table 1 :

 1 Proximal operators of γf for different nonlinear functions f (see Sect. 2.5.2)

	RR n°8837

  d is the adjoint operator of × d (see Appendix B for definition and properties). From Lemma B.1 and Lemma B.2, we have s× d . 2 ≤ s 1 . One can easily verify that in 3D, H 2 ≤ 12. Hence the Lipschitz constant of ∇g 1 is bounded by

RR n°8837

† † Let us emphasize again that estimating receptive fields is an important step done in almost every recording session. However, it is generally not the main goal of the experimentalist who want to investigate other properties of cells. Thus, because of physiological constraints such as the limited life duration of the tissue, keeping the receptive field estimation step short will let more time to experimentalists for the main purpose of their main RR n°8837

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements We thank Daniela Pamplona (Inria, Biovision team) as well as Evelyne Sernagor and Gerrit Hilgen (Institute of Neuroscience, Newcastle University, UK) and Matthias Hennig (Institute for Adaptive and Neural Computation, University of Edinburgh) for their helpful insights.

Appendices

A Bayesian approach

This section justifies the formulation (3) derived from (P) using a Bayesian approach. Before giving the proof, some definitions and properties are reminded.

Definition A.1. Let λ(t) ∈ C 0 (R + , R + ). A process that produces random points in time is a Poisson random process of rate λ(t) if the counting process n(t) (number of events at t) satisfies the following properties ∀t ≥ 0 :

The arrival time t (i) for i ≥ 1 are defined by

Lemma A.1. Let n(t) be a Poisson random process of rate λ(t), t 0 > 0 and t ≥ t 0 , the probability that n events have occurred from time t 0 to time t is

where

λ(s)ds.

study.

Inria

E.1 Computation of y 1 (x)

.

The equation has two solutions which are:

and ϕ x,γ varies as

We search values of x for which y ± are in [- 1 2 , 1 2 ]:

Since γ ≤ 1 6 , we easily check that y + ≥ 1 2 for all x and that x 0 ≤ 1 2 . For y -there are three cases:

and so y 1 (x) is given by

, where y -is given by [START_REF] Wilkie | Model completeness results for expansions of the ordered field of real numbers by restricted Pfaffian functions and the exponential function[END_REF].

Case 2: ∆ < 0 (i.e. x > x 0 ) In this case, ϕ x,γ ≤ 0 i.e. ϕ x,γ is non increasing and y 1 (x) = 1 2 .

E.2 Computation of y 2 (x)

By a similar reasoning as for the computation of g 1 and y 1 , we get

and from ( 43), we deduce that ϕ x,γ (y) is minimal at y = y 1 (x) and its value is g 1 (x). Hence g 1 (x) ≤ min(ϕ x,γ (- 1 2 ), ϕ x,γ ( 1 2 )) = g 2 (x) and g(x) = g 1 (x) and prox γf (x) = y 1 (x).

(ii) If x ≥ 1 2 and x ≤ -1 2 it is easily seen that g(x) = g 2 (x) and prox γf (x) = y 2 (x) = x.

To summarize, prox γf (x) is given by (Fig. 14): Remark E.2. We remind that the case γ ≥ 1 6 never happens in our applications since γ is small since it is a O 1 α and α is large. However in this case we can show that the proximal operator is multivalued. More precisely, there exists x 1 ∈ [ 1 2 , x 0 ] such that y -(x 1 ) < x 1 and the proximal operator is for γ ≥ 1 6 of the form:

and so prox γf (x 1 ) = {y -(x 1 ), x 1 } i.e. the proximal operator is multivalued.