Sylvain Rubenthaler

OPTIMAL HEDGING IN DISCRETE TIME

, we present discrete time formulas minimizing the mean square hedging error for multidimensional assets. In particular, we give explicit formulas when a regime-switching random walk or a GARCH-type process is utilized to model the returns. Monte Carlo simulations are used to compare the optimal and delta hedging methods.

Hedging, option pricing, GARCH, regime-switching

Introduction

In many applications, one is interested in finding a discrete-time dynamically traded portfolio so that its value at maturity is as close as possible to a target function of the underlying assets. When the target function is a payoff, this can be interpreted as option pricing and hedging. However, sometimes the target function is not a payoff as it happens when one tries to replicate hedge funds or create synthetic funds [START_REF] Papageorgiou | Replicating the properties of hedge fund returns[END_REF]. Assuming that the error measure is the average quadratic hedging error, [START_REF] Schweizer | Variance-optimal hedging in discrete time[END_REF] solved the hedging problem for one risky asset. He showed that the initial value of the portfolio, which can be interpreted as the "value" of the option, is the average, under the "real probability measure", of the discounted payoff, multiplied by a martingale, which is not necessarily positive. In the latter case, the martingale cannot be used as the density of an equivalent martingale measure.

Partial funding in support of this work was provided by the Natural Sciences and Engineering Research Council of Canada, by the Fonds Québécois de Recherche sur la Nature et les Technologies, the Institut de Finance Mathématique de Montréal and by the PPF Complexité-Modélisation-Finance de l'Universtité Nice-Sophia Antipolis. We would like to thank Gerasimos Rassias for his helpful suggestions.

Even if the hedging problem has been solved quite generally by [START_REF] Schweizer | Variance-optimal hedging in discrete time[END_REF] in the one-dimensional case, it seems to have been ignored or forgotten, e.g., [START_REF] Bouchaud | Back to basics: historical option pricing revisited[END_REF] or [START_REF] Cornalba | Option pricing and hedging with temporal correlations[END_REF]. More troubling, delta hedging, based on the Black-Scholes-Merton model, is still used in practice even if it has been shown that the geometric Brownian motion model is an inadequate model for the underlying assets [START_REF] Kat | Who needs hedge funds? A copula-based approach to hedge fund return replication[END_REF]. Even when the geometric Brownian motion model is adequate, the hedging error in discrete time is not zero, converging only to zero as the number of hedging periods tends to infinity. See, e.g., [START_REF] Boyle | Discretely adjusted option hedges[END_REF], Wilmott (2006, Chapters 46-47). In addition, when the market is not complete but there is no arbitrage, there are infinitely many equivalent martingale measures (EMM). One then has to choose the "best" martingale measure with respect to some utility criterion. There is a huge literature on this subject. Indeed, when a NGARCH process is utilized to model the log-returns, [START_REF] Duan | The GARCH option pricing model[END_REF] proposes a solution to the EMM problem. Unfortunately, Duan also suggests a delta hedging strategy, which has been shown to be wrong by [START_REF] Garcia | A note on hedging in ARCH and stochastic volatility option pricing models[END_REF].

Motivated by applications in hedge fund replication, [START_REF] Papageorgiou | Replicating the properties of hedge fund returns[END_REF] proposed a locally optimal solution minimizing the average quadratic hedging error at each period for the general multidimensional asset case. They erroneously claimed that it was globally optimal, which is only true if the discounted underlying assets are martingales. A first motivation for the present paper is to correct that mistake and give explicit formulas for the results in [START_REF] Černý | On the structure of general mean-variance hedging strategies[END_REF] and generalizing those of [START_REF] Schweizer | Variance-optimal hedging in discrete time[END_REF]. A second motivation is to show that, when regime-switching random walks and GARCH processes are used to model the returns, the optimal solution of the hedging problem yields superior results to those obtained by delta hedging.

The optimal solution of the discrete time hedging problem is described in Section 2, giving explicit expressions for the results of [START_REF] Černý | On the structure of general mean-variance hedging strategies[END_REF]. It is worth noting that when the asset value process is Markovian, or a component of a Markov process, the optimal solution can be implemented using approximation techniques of dynamic programming. Two such cases are considered. Finally, in Section 3, simulations are used to compare optimal hedging with delta hedging for geometric random walks and NGARCH models.

Optimal hedging strategy in discrete time

Denote the price process by S, i.e., S k is the value of the d underlying assets at period k and let F = {F k , k = 0, . . . , n} be a filtration under which S is adapted.

Assume that S is square integrable. Set ∆ k = β k S k -β k-1 S k-1 , where the discounting factors β k are predictable, i.e. β k is F k-1 -measurable for k = 1, . . . , n.
The aim of this section is to find an initial investment amount V 0 and a pre-

dictable investment strategy -→ φ = (φ k ) n k=1 such that φ ⊤ k ∆ k is square integrable and which minimizes the expected quadratic hedging error E G V 0 , - → φ 2
, where

G = G V 0 , - → φ = β n C -V n , and V k = V 0 + k j=1 φ ⊤ j ∆ j , k = 0, . . . , n.
The existence of an optimal solution was proven first in the univariate case by [START_REF] Schweizer | Variance-optimal hedging in discrete time[END_REF]. [START_REF] Motoczyński | Multidimensional variance-optimal hedging in discrete-time model-a general approach[END_REF] considered a multivariate setting without furnishing explicit solutions. Finally, [START_REF] Černý | On the structure of general mean-variance hedging strategies[END_REF] treated a much more general case. However, in the discrete time case, it is faster and easier to find directly the explicit formulas than trying to recover them from their results.

2.1. Offline computations. Once a dynamic model is chosen for the asset prices, one must start with some computations that are necessary for the implementation.

Set P n+1 = 1, γ n+1 = 1, and for k = n, . . . , 1, define

A k = E ∆ k ∆ ⊤ k P k+1 |F k-1 , µ k = E (∆ k P k+1 |F k-1 ), b k = A -1 k µ k , P k = n j=k 1 -b ⊤ j ∆ j , and γ k = E(P k |F k-1
), provided these expressions exist. Under some extra assumptions given below, it can be shown that they are indeed well defined. The proof is given in Appendix A.1.

Lemma 2.1.1. Suppose that E(γ k+1 |F k-1 )A k -µ k µ ⊤
k is invertible P-a.s., for every k = n, . . . , 1. Then γ k ∈ (0, 1] and A k is invertible for all k = 1, . . . , n. In addition

(γ k+1 ) n k=0 is a positive submartingale.
Remark 2.1.2. In the univariate case, [START_REF] Schweizer | Variance-optimal hedging in discrete time[END_REF] states sufficient conditions for the validity of the assumptions of Lemma 2.1.1. It is not obvious how they could be generalized to the multivariate case. Therefore, in most applications, one has to verify these conditions, often using brute force calculations.

2.2.

Optimal solution of the hedging problem.

Theorem 2.2.1. Under the assumptions of Lemma 2.1.1, the solution

V 0 , - → φ of the minimization problem is V 0 = E(β n CP 1 )/γ 1 , and φ k = α k -V k-1 b k , where α k = A -1 k E (β n C∆ k P k+1 |F k-1 ), k = n, . . . , 1.
The proof is given in Appendix A.2.

2.2.1. Option value. Let C k be the optimal investment at period k, so that the value of the portfolio at period n is as close as possible to C, in terms of mean square error. One could then interpret C k as the value of the option at period k, for any k = 0, . . . , n. It then follows from Theorem 2.2.1 that C k is given by (2.2.1)

β k C k = E(β n CP k+1 |F k ) E(P k+1 |F k ) , k = 0, . . . , n, so one can write (2.2.2) β k-1 C k-1 = 1 γ k E β k C k 1 -b ⊤ k ∆ k γ k+1 |F k-1 = E (β n CU k • • • U n |F k-1 ) ,
where

U k = E(P k |F k ) E(P k |F k-1 ) , k = 1, . . . , n + 1, while an alternative expression for α k is (2.2.3) α k = A -1 k E (β k C k ∆ k γ k+1 |F k-1 ) . Remark 2.2.2. Setting Z 0 = 1 and Z k = k j=1 U j , k = 1, . . . , n, one obtains that (Z k , β k C k Z k , β k S k Z k ) n
k=0 are martingales. However, in most applications, Z does not define a change of measure unless it takes only positive values.

2.2.2. Implementation issues. If the process S is Markov and

C n = C n (S n ), then C k = C k (S k ), α k = α k (S k-1 ), and b k = b k (S k-1
). It follows that all these functions can be approximated using the methodology developed in [START_REF] Papageorgiou | Replicating the properties of hedge fund returns[END_REF].

Another interesting case is when S k is not a Markov process but (S k , h k ) is, even if h k is not observable, as in GARCH and regime-switching models. In this case, if

C n = C n (S n ), then C k = C k (S k , h k ), α k = α k (S k-1 , h k-1 ), b k = b k (S k-1 , h k-1
), and .

γ k = γ k (S k-1 , h k-1
), for k = 1, . . . , n + 1. All these functions can be approximated using the methodology developed in [START_REF] Rémillard | Option Pricing and Dynamic Discrete Time Hedging for Regime-Switching Geometric Random Walks Models[END_REF] for the regime-switching case. Implementation of the hedging strategy then requires predicting h t .

Remark 2.2.3. One could suggest to use the smallest filtration to get rid of the unobservable process h but in this case, all conditional expectations based on F k would depend on all past values S 0 , . . . , S k , making it impossible to implement in practice.

Verification of the assumptions of Lemma 2.1.1.

In what follows, we consider some interesting models used in practice, for which it is possible to show that the assumptions of Lemma 2.1.1 hold true and that the optimal solution can be computed via a dynamic program.

2.3.1. Regime-switching geometric random walks. An interesting model, which includes geometric random walk models, is to consider a regime-switching geometric random walk. Theses models can display serial dependence in the log-returns and may account for changing volatility over time. For implementation issues, including estimation and goodness-of-fit tests, see, e.g., [START_REF] Rémillard | Option Pricing and Dynamic Discrete Time Hedging for Regime-Switching Geometric Random Walks Models[END_REF].

To define the process, suppose that τ is a finite homogeneous Markov chain with transition matrix Q with values in {1, . . . , l} representing the non-observable regimes and set

β k S i k = S i 0 k t=1 {1 + ξ i t }, i = 1, . . . , d
, where, given τ 1 = i 1 , . . . , τ n = i n , ξ 1 , . . . , ξ n are independent with ξ j ∼ P i j , j = 1, . . . , n, E i (ξ j ) = E(ξ j |τ j = i) = m(i),

and E i ξ j ξ ⊤ j = B(i).
The interpretation of the model is easy: At a given period t, a regime τ t is chosen at random, according to the Markov chain model, and given τ t = i, ξ t is chosen at random according to distribution P i . When there is only one regime, one obtains a geometric random walk where all ξs are independent.

We assume that the B(i) -m(i)m(i) ⊤ is invertible for any i = 1, . . . , l. Setting

X k = β k S k , one gets ∆ k = X k -X k-1 = D(X k-1 )ξ k , k = 1, . . . ,
= γ k (τ k-1 ), γ k (i) ∈ (0, 1], A k = A k (S k-1 , τ k-1 ) and b k = b k (S k-1 , τ k-1 ), where (2.3.1) A k (s, i) = β 2 k-1 D(s) l j=1 Q ij γ k+1 (j)B(j) D(s), (2.3.2) b k (s, i) = D -1 (β k-1 s)ρ k+1 (i), (2.3.3) γ k (i) = l j=1 Q ij γ k+1 (j) 1 -ρ k+1 (i) ⊤ m(j) , with ρ k+1 (i) = l j=1 Q ij γ k+1 (j)B(j) -1 l j=1 Q ij γ k+1 (j)m(j) . If in addition C = C(S n ), then C k = C k (S k , τ k ) and α k = α k (S k-1 , τ k-1
), where

C k-1 (s, i) = β k β k-1 l j=1 Q ij γ k+1 (j) γ k (i) × C k β k-1 β k D(s)(1 + y), j 1 -ρ k+1 (i) ⊤ y P j (dy), α k (s, i) = β k β k-1 D -1 (s) l j=1 Q ij γ k+1 (j)B(j) -1 l j=1 Q ij γ k+1 (j) × C k β k-1 β k D(s)(1 + y), j yP j (dy).
2.3.2. GARCH-type models. Here, one assumes that

∆ k = β k S k -β k-S k-1 = β k-1 S k-1 ξ k , with ξ k = π 1 (h k-1 , ǫ k ), and h k = π 2 (h k-1 , ǫ k )
= γ k (h k-1 ) and A k (s, h) = β 2 k-1 s 2 B k (h), b k (s, h) = m k (h) sβ k-1 B k (h) , where B k (h) = π 2 1 (h, y)γ k+1 {π 2 (h, y)} ν(dy), m k (h) = π 1 (h, y)γ k+1 {π 2 (h, y)} ν(dy), and γ k (h) = 1 -m k (h) B k (h) π 1 (h, y) γ k+1 {π 2 (h, y)} ν(dy). Also, if C = C n (S n ), then C k-1 (s, h) = β 1 γ k (h) C k s β 1 {1 + π 1 (h, y)}, π 2 (h, y) ×γ k+1 {π 2 (h, y)} 1 - m k (h) B k (h) π 1 (h, y) ν(dy), α k (s, h) = β 1 sB k (h) C k s β 1 {1 + π 1 (h, y)}, π 2 (h, y) ×γ k+1 {π 2 (h, y)} π 1 (h, y)ν(dy).
Hence, the optimal solution can be written as a dynamic program.

Examples of application

In this section we consider pricing and hedging of European calls for two geometric random walk models, when the returns are i.i.d. Gaussian and i.i.d. differences of Laplace distributions, and for a NGARCH model. It follows from the previous sections that optimal hedging solutions exist for these cases, and the optimal solution can be written as a dynamic program associated with functions of a finite number of variables.

For solving such dynamic programs, we discretize the state space into a finite grid and we compute approximations of expectations using Monte Carlo simulations at every point of the grid. Linear interpolations are used for points outside the grid at each time step. Since the expectations are always with respect to the same probability measure, only one sequence of random numbers may be used, using the ideas in Del [START_REF] Del Moral | Monte Carlo approximations of American options[END_REF][START_REF] Del Moral | Monte Carlo Approximations of American Options that Preserve Monotonicity and Convexity[END_REF].

3.1. Geometric random walk models. Here we consider discretized versions of the Black-Scholes (BS) and Variance Gamma (VG) models for the underlying asset over 23 periods. In each case, the 22 periodic returns are i.i.d., so that the mean and volatility at maturity are respectively 9% and 6%. For the BS model, the returns are Gaussian, while for the VG model, the returns are differences of i.i.d. Gamma variates, so that the distribution at maturity is Laplace (double exponential). These models are particular cases of regime-switching models with only one regime. We do not consider regime-switching models since it has been done in [START_REF] Rémillard | Option Pricing and Dynamic Discrete Time Hedging for Regime-Switching Geometric Random Walks Models[END_REF], where the daily log-returns of the S&P 500 are analyzed.

We are going to price and hedge a call with strike K = 100 and maturity 1 year, using 22 replication periods and a 2000 points discretization of the asset values over the interval [80,120]. The annual rate is 5%. A sequence of 10000 random points were used for the computation of functions α k and A k , while 10000 paths were used to compute the hedging errors. Delta hedging is optimal in the continuous time limit for the BS model, but not for the VG model. As expected, according to Figures 12, the values of the call C 0 and initial investment strategy φ 1 , obtained from the optimal hedging, are close to those obtained using the Black-Scholes formula (even with 22 hedging periods), while they differ for the VG model. This is also reflected in the distribution of the hedging errors, as illustrated in Figure 3. Descriptive statistics of the hedging errors are given in Table 1. Simulations can also be used to show that as the number of hedging periods increases, the hedging error tends to zero for the BS model, while it is never 0 for the VG model. Note that the RMSE of the optimal hedging is always less than the one of the delta hedging.

ξ k -1 = r + λ h k-1 -1 2 h k-1 + h k-1 ε k , and h k = α 0 + α 1 h k-1 ε 2 k + β 1 h k-1
, with ε k ∼ N(0, 1) and parameters α 0 = 1.524 × 10 -5 , α 1 = 0.1883, β 1 = 0.7162 and λ = 7.452 × 10 -3 . Under the EMM, we have e

ξ k -1 = r + -1 2 h k-1 + h k-1 ε k and h k = α 0 + α 1 h k-1 (ε k -λ) 2 + β 1 h k-1 .
We price and hedge a call with strike K = 100 and maturity 30 days using daily replication, using a grid of 500 points for the asset on [60,140], while the grid for the volatility consists in 90 points of the interval [.00005, .0007]. The annual rate is 0%. In what follows, B&S hedging means delta hedging using the B&S formulas, while Duan's methods consists in picking his suggested EMM and taking the delta of the option. The value of the option and the initial number of asset are displayed in Figure 4, while descriptive statistics of the 10000 hedging errors are given in Table 2 for the three hedging methodologies, showing that the errors are more concentrated about 0 for the optimal hedging.

Conclusion

In this paper we presented the optimal solution for a discrete time hedging portfolio. When the underlying process is Markov or a component of a Markov process, the optimal hedging strategy depends on deterministic functions that can be approximated. We also find explicit formulas for two interesting models. Finally, numerical simulations show that optimal hedging is preferable to delta hedging.

µ ⊤ n A -1 n µ n = µ ⊤ n Σ -1 n µn 1+µ ⊤ n Σ -1
n µn , so γ n = A.2. Proof of Theorem 2.2.1. Using the proof of Lemma 2.1.1, one can easily check that a k , b k and φ k make sense and that φ ⊤ k ∆ k is square integrable. Next, it is easy to check that a necessary and sufficient condition for V 0 , -→ φ to minimize

E G V 0 , -→ φ 2 is that E G V 0 , -→ φ = 0 and E G V 0 , - → φ ∆ k |F k-1 = 0 for
all k = 1, . . . , n. The necessity comes from the fact that for any event O ∈ F k-1 , one must have 0 = d dǫ ǫ=0 E G V 0 ,

- → φ -ǫI O ∆ k 2 = -2E G V 0 , - → φ ∆ k I O ,
which is equivalent to the condition E G V 0 , -→ φ ∆ k |F k-1 = 0, while the condition E G V 0 , -→ φ = 0 comes from the fact that for any θ, one must have

0 = d dǫ ǫ=0 E G V 0 + ǫθ, -→ φ 2 = -2E G V 0 , - → φ .
To see that the conditions are sufficient, it suffices to check that

E G V 0 + θ 0 , ---→ φ + ψ 2 = E G V 0 , - → φ 2 + E    θ 0 + n k=1 ψ ⊤ k ∆ k 2    .
The proof that -→ φ is the solution is based on the following equation, which can be easily proven by induction.

(A.2.1)

E(V n |F k ) = V k E(P k+1 |F k ) + E {β n C(1 -P k+1 )|F k } , k = 1, . . . , n.
To complete the proof of theorem, note that from (A.

  n, where D(s) is the diagonal matrix constructed from vector s. Note that S is not a Markov process in general but (S, τ ) is a Markov process. The validity of the assumptions of Lemma 2.1.1 follows from the next result, proved in Appendix A.3. Proposition 2.3.1. For any k = 1, . . . , n and i = 1, . . . , l, γ k

Figure 1 .

 1 Figure 1. Call option value C 0 for the Black-Scholes (left panel) and Variance Gamma models (right panel) with 22 periods of hedging.

Figure 2 .

 2 Figure 2. Initial investment strategy φ 1 in the underlying asset for the Black-Scholes (left panel) and Variance Gamma models (right panel) with 22 periods of hedging.

Figure 3 .

 3 Figure 3. Estimated densities of the hedging error G for the Black-Scholes (left panel) and Variance Gamma models (right panel) with 22 periods of hedging.

Figure 4 .

 4 Figure 4. Optimal hedging call option value C 0 and initial investment strategy φ 1 for the NGARCH model with 30 periods.

n

  µn ∈ (0, 1]. As a result, γ n ≤ E(γ n+1 F n ) = 1.The rest of the proof follows easily by reverse induction, using Proposition B.0.1 with the mean and covariance matrix of ∆ k under the probability distribution Q k , withQ k (O) = E(I O γ k+1 |F k-1 )/E(γ k+1 |F k-1 ), O ∈ F k , for k = n -1, . . . , 1.

  2.1), (A.2.2) E G V 0 , -→ φ |F k = E(β n CP k+1 |F k ) -V k E(P k+1 |F k ), k = 0, . . . , n.

Table 1 .

 1 Statistics of hedging errors (Payoff-Portfolio) for the Black-Scholes and Variance Gamma models.

		Black-Scoles	Variance Gamma
	Stats	Optimal	Delta Optimal	Delta
	Average	-0.0065 0.0076	0.0151 0.0518
	Median	-0.0014 0.0029 -0.1543 -0.2886
	Volatility	0.2537 0.2774	1.0510 1.2529
	Skewness	1.3781 0.5515	6.1479 6.0965
	Kurtosis	22.0975 8.8391 63.3897 60.9747
	Minimum	-1.7986 -1.7578 -4.4372 -1.5160
	Maximum	3.5389 2.2978 18.5163 23.6422
	VaR(99%)	0.8231 0.9244	4.6976 5.9450
	VaR(99.9%)	1.9576 1.5779 11.7905 12.2222
	RMSE	0.2538 0.2775	1.0511 1.2540
	3.2. NGARCH model. As in Duan (1995), we consider the NGARCH model
	where e

Table 2 .

 2 Statistics of hedging errors (Payoff-Portfolio) for the NGARCH model.

	Stats	Optimal	Delta	Duan
	Average	-0.0159 -0.0954 0.0085
	Median	-0.1549 -0.2094 -0.1297
	Volatility	0.8568 0.8951 0.9032
	Skewness	1.7205 2.7558 2.9171
	Kurtosis	10.5790 29.3709 29.9947
	Minimum	-1.9966 -2.2302 -2.0789
	Maximum	9.9114 17.4873 17.7232
	VaR(99%)	2.7698 2.8613 3.0831
	VaR(99.9%)	5.4893 6.3574 6.5070
	RMSE	0.8569 0.9001 0.9033

Appendix A. Proofs of the main results

A.1. Proof of Lemma 2.1.1. First, we will show that γ k ∈ (0, 1] and A k is invertible for all k = 1, . . . , n. By hypothesis, E(γ k+1 |F k-1 )A k -µ k µ ⊤ k is invertible for all k = 1, . . . , n. In particular, it is true for k = n, yielding that Σ n = A n -µ n µ ⊤ n is invertible, which is the conditional covariance matrix of ∆ n given F n-1 . It then follows from Proposition B.0.1 that A n is invertible. Without loss of generality, one may assume that A n is diagonal. Otherwise, we diagonalize it in the form

A.3. Proof of Proposition 2.3.1. The result is obviously true for k = n + 1.

Suppose it is true for k + 1. For i given, set

. By hypothesis, π 1 , . . . , π l are probabilities adding to 1, so if X ∼ P j with probability π j , then γ k (i) = D 1 -µ ⊤ B -1 µ , where µ = E(X) and B = E XX ⊤ . Let Σ be the covariance matrix of X. It is non singular since the covariance of X under P j is non singular. It then follows from Proposition B.0.1 that

Since D > 0 by hypothesis, one may conclude that γ k (i) > 0. As a by-product we get that γ k (i) ≤ 1 if γ k+1 (j) ≤ 1 for all j = 1, . . .. Since that is true for γ n+1 ≡ 1, one may conclude that for all k = 1, . . . , n, γ k (i) ≤ 1.

The rest of the proof is easy.

Appendix B. Auxiliary results

Proposition B.0.1. Suppose A = Σ+bb ⊤ where Σ is symmetric and invertible. Then A is invertible, and

Service de l'enseignement des méthodes quantitatives de gestion and GERAD, HEC Montréal, 3000 chemin de la Côte Sainte-Catherine, Montréal (Québec) Canada

H3T 2A7

E-mail address: bruno.remillard@hec.ca Laboratoire J.-A. Dieudonné, Université de Nice-Sophia Antipolis, Parc Valrose, 06108 Nice Cedex 2, France E-mail address: rubentha@unice.fr


Bruno Cessac

Rodrigo Salas

Thierry Viéville

Using event-based metric for event-based neural network weight adjustment

The problem of adjusting the parameters of an event-based network model is addressed here at the programmatic level. Considering temporal processing, the goal is to adjust the network units weights so that the outcoming events correspond to what is desired. The present work proposes, in the deterministic and discrete case, a way to adapt usual alignment metrics in order to derive suitable adjustment rules. At the numerical level, the stability and unbiasness of the method is verified.

Introduction

Studying the computational power of neural networks with event-based activity (e.g.: [START_REF] Rieke | Spikes, Exploring the Neural Code[END_REF][START_REF] Gerstner | Spiking Neuron Models[END_REF]) is a well-addressed topic, see [START_REF] Paugam-Moisy | Handbook of Natural Computing, chapter Computing with Spiking Neuron Networks[END_REF][START_REF] Schrauwen | Towards Applicable Spiking Neural Networks[END_REF] for a recent review about spiking network computation, while [START_REF] Bruno Cessac | Overview of facts and issues about neural coding by spikes[END_REF] provides a detailed discussion on temporal aspects of such computations. See [START_REF] Viéville | Using event-based metric for network weight adjustment (extended version)[END_REF] for further details on the related modeling choices. In order to contribute to this general topic, we develop here a framework allowing us to effectively adjust the network parameters in order to tune the outcoming events.

Position of the problem

We consider an input/output dynamical system with N units, governed by a recurrent function, V = {• • • V n [t] • • • } being the output state variable value of the units of output index n ∈ {0, N { (i.e., 0 ≤ n < N ) at time t ∈ {0, T {. Some output units may be "hidden", i.e. not observed. Here W stands for the network parameters or "weights", to be tuned. The exact form of V is not relevant at this stage, but the gradient ∇ W V W must be welldefined in order to adjust W. One track is to consider regular forms of V. For a spiking neuron network this means that we have to consider either Hodgkin-Huxley equations, or some suitable reduction like the FitzHugh-Nagumo model or the SRM model [START_REF] Gerstner | Spiking Neuron Models[END_REF]. Another track, is to "mollify" V, i.e., defines it as the limit of a series of regular functions, as experimented in [START_REF] Viéville | A general algorithm to estimate recurrent networks weights (extended version)[END_REF].

We define an event

Z n [t] def = H(V n [t] -θ) ∈ {0, 1}
, where H is the Heaviside function, as the fact that the output value is higher than a threshold θ. The goal is thus to adjust the output events Z of the deterministic discrete-time dynamical system V, with respect to a reference output events Z. The key point is to deal with the fact that the notion of event is intrinsically "discontinuous".

Considering alignment metric We define the distance between two finite eventtrains Z, Z as the minimum cost of transforming one event-train into another. See [START_REF] Victor | Spike train metrics[END_REF][START_REF] Houghton | Measuring representational distances -the spike-train metrics approach[END_REF] for a general introduction. Following [START_REF] Bruno Cessac | Overview of facts and issues about neural coding by spikes[END_REF], we consider a generalized alignment metric: non-stationary cost (e.g., recent events may count more than older ones) and non-linear shift (e.g., neglecting tiny delays), as described in Fig. 1. Two kinds of operations are defined for an alignment metric. (i) Event insertion/deletion, the cost of each operation being set to γ ± t at time t, e.g., γ ± t = 1, while non-stationary different insertion/deletion costs may be defined.

(ii) Event shift, the cost to shift from one event in Z at time t to one event Z at time t, being an increasing positive function of the non-stationary normalized shift delay φt(( t-t)/τ ), for a given time-constant τ (e.g. φt(( t-t)/τ ) = | t-t|/τ ), while non-stationary non-linear different forward/backward shift-cost may be defined, since φ() is parameterized by t.

From the upper to the lower event train is shown, using from top to bottom an insertion, a rightward shift, a leftward shift and a deletion respectively Fig. 1: An example of minimal alignment (borrowed from [START_REF] Victor | Spike train metrics[END_REF]).

Obviously the distance is zero (no editing operation) if and only if both trains are equal, is always bounded by the number of events in both event-trains (i.e. the cost of deleting/inserting all events), thus also by twice the number of samples in the discretized case. For small τ , the distance approaches the number of non-coincident events, since instead of shifting events it is cheaper to insert/delete non-coincident events, while when τ → 0, γ ± t = 1 we obtain the coincidence (or Hamming) distance equal to the number of non-coincident events. Given two time sequences with the same number of events, there is always a τ high enough for the distance to correspond to the weighted sum of time differences between both train events, as used in, e.g., [START_REF] Schrauwen | Towards Applicable Spiking Neural Networks[END_REF]. More generally, for high τ , the distance basically equals the difference in event number (rate distance) [START_REF] Victor | Spike train metrics[END_REF]. When considering event-trains with more than one unit, our approach consists to sum the distances for each unit alignment, i.e., consider each unit independently, avoiding the related estimation to suffer from NP-completeness [START_REF] Houghton | Measuring representational distances -the spike-train metrics approach[END_REF].

Since we want to tune the Z events in order to approximate the Z events, let us introduce an alignment indexation as follow: δ(t) = t-t if the two events Z[t] = 1 and Z[ t] = 1 are aligned by a shift, δ(t) = ±0 to code for an insertion/deletion, while δ(t) = 0 otherwise. In words, we not only compute the distance but make explicit the alignment operations (shif t, deletion, insertion) allowing to "edit" Z in order to obtain Z. The δ code function is used to explicitly match both trains. The distance dk ,k between the first k events in Z and the first k events Z and the related δ indexing are iteratively defined by induction, after [START_REF] Victor | Spike train metrics[END_REF] but now generalized (see [START_REF] Viéville | A general algorithm to estimate recurrent networks weights (extended version)[END_REF] for a detailed derivation). We write t k , k > 0 the k-th value such that Z[t k ] = 1, with a similar notation for tk, k > 0. On one hand, dk ,0 = l< k γ - tl , since the distance between any event-train and the empty event-train corresponds to the cost of deleting all events, while δ(tk) = -0 in this case. Similarly, d 0,k = l<k γ + t l corresponds to inserting all events, with δ(tk) = +0. On the other hand:

dk +1,k+1 = min      dk ,k+1 + γ - tk (deletion), δ(t k+1 ) = -0 dk +1,k + γ + t k (insertion), δ(t k ) = +0 dk ,k + φt k tk-tk τ (shif t), δ(t k ) = tk -t k . (1)
Obviously, several alignment operation sequences may lead to the same minimal alignment cost. In order to make a choice, from the last time to the previous time, we consider that shift is preferable to insertion/deletion, since it is a reasonable assumption to heuristic that it is going to have a less important influence on the dynamics than the apparition/cancellation of an unexpected event. This defines algorithmically a unique well-defined indexing function for a given distance, as illustrated in Fig. 2, thus solves the ambiguities. On the reverse, solving these ambiguities allows us to define algorithmically a unique indexing function.

Although computing such a distance and indexes seems subject to a combinatorial complexity, this is a quadratic algorithm (i.e. with a complexity equal to the product of the numbers of events), and its derivation, done by induction, is similar to usual alignment distance calculations [START_REF] Houghton | Measuring representational distances -the spike-train metrics approach[END_REF][START_REF] Bruno Cessac | Overview of facts and issues about neural coding by spikes[END_REF]. Regarding indexing, this means that we do not have to explore all possible alignment operation sequences, in order to define a globally well-defined process. This indexing definition also allows us to enrich the original alignment distance by not only considering a composite number describing the distance in terms of shift and insertion/deletion, but allowing to make explicit a numerical approximation of the number of shifts versus insertion/deletion. This is the same feature as in message-passing alignment mechanisms [START_REF] Dauwels | Measuring neural synchrony by message passing[END_REF], but defined here in a much simpler context.

Edition is realized from the most recent event to the oldest event in the past, i.e., from right to left. The previous specification solves the ambiguity. A From right to left, shift is preferred to insertion, thus shift precedes insertion. B From right to left, two minimal shifts are preferred to a higher shift and a event coincidence.

Mollified version of the alignment distance

From the previous construction we now introduce the key idea of the paper, i.e., propose a variational expression of the alignment distance. To this aim, we "mollify" the event generation mechanism, i.e. replace the Heaviside function by a suitable convolution

H υ = υ * H = H(u+ √ υ ν) exp(-υ/(u+ √ υ ν))
where ν is a margin maintaining the state at a non-infinitesimal distance to the threshold, and υ → 0 is the mollification factor (see [START_REF] Viéville | Using event-based metric for network weight adjustment (extended version)[END_REF] for details)), as show in Fig. 3. We obtain after some algebra given in [START_REF] Viéville | Using event-based metric for network weight adjustment (extended version)[END_REF]: d Z, Z = lim υ→0 d υ Z, V , with:

d υ Z, V = nt,δ[t]=0 γ ± t H υ (1 -2 Zn [t]) (V n [t] -θ) + nt,δ[t] =0 φ t δ[t] τ Z n [t]+ (1 -Z n [t]) H υ -δ[t] τ H υ (θ -V n [t]) ν=0 (2)
The key point is that now the criterion is not defined with respect to Z but V. Qualitatively an increase of V n [t] tends to shift event in the past, avoid deletion but induce insertion of event, whereas a decrease of V n [t] tends to shift event in the future, induce deletion but avoid insertion of event. Changes are now differentiable, thanks to the mollification and the gradient

∇ W d υ Z, V = ∇ V d υ Z, V ∇ W V
is obvious to derive. We thus can tune V, thus W to optimize the alignment metric as desired, with a straight-forward implementation for a feed-forward system and the need of specific method for a recurrent structure, as developed elsewhere [START_REF] Viéville | A general algorithm to estimate recurrent networks weights (extended version)[END_REF]. The metric allows us to calculate the appropriate network weights by a simple numerical minimization. Since H υ is convex for suitable υ [START_REF] Viéville | Using event-based metric for network weight adjustment (extended version)[END_REF], for a fixed value of δ the criterion is convex as the sum of positively weighted functions. However the criterion is also optimized with respect to δ and as soon as an event occurrence is modified by a variation of the weights, the indexing is to be recalculated, while, up to our best knowledge, there is no chance to guarantee a global minimum, so that we now turn to numerical verification.

It is drawn here for ν = 0 and in black, brown, red, orange, yellow, green, blue, for υ = [1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01], respectively. The curves are convex below the magenta horizontal line.

Numerical experiment

Criterion calibration. In order to estimate the performance of the estimation, we must first quantify to which extends we obtain numbers "better than by chance", i.e., to which extends the minimized alignment distance yields a better result than if the output would have been random (a small distance may simply mean that events are sparse !). In order to obtain a correct order of magnitude, we considered the normalized alignment distance with respect to a random event train of the same rate. More precisely, if two event trains of T samples are drawn from a Bernoulli distribution of rates r (i.e., samples are random and independent), it is straightforward to obtain the average coincidence distance, i.e. E[d( Z, Z)] = 2 T r (1-r). However, the same derivation is not obvious for an alignment distance parameterized with τ and we have numerically interpolated the value as a power of τ , for the standard alignment distance, obtaining for γ ± t = 1, φ t (s) = |s|, with a residual standard deviation better than 1.5%:

E[d( Z, Z)] = 2 T r (1 -r) 1.183 + 0.183 r (1 -r) (τ + 1) 0.265+1.444 r (1-r) .
Numerical robustness. Let us now illustrate the previous developments considering a leaky integrate and fire (LIF) network, as in [START_REF] Rostro-Gonzalez | Backengineering of spiking neural networks parameters[END_REF]. As a test, we have generated hundred of input/output data sets using a "master" network and have verified that the learning algorithm applied on another model of the same dimensions is able to find weights that reproduce the input/output function.

Weights values are randomly drawn from a Gaussian distribution of zero mean and standard-deviation σ ∈ [0.1, 10]. The LIF resetting mechanism is mollified as for the event thresholding. This is a basic verification of both the correctness of the code and the numerical stability of the estimation. Hundred of runs have provided correct results, as expected. For long length T > 10 5 and complex dynamics the method may fail finding the exact solution with the standard parameters. For small length epoch, as expected, there is always an exact solution, in fact there is one, even if the raster is not generated from the same model [START_REF] Rostro-Gonzalez | Backengineering of spiking neural networks parameters[END_REF]. These tests have been done for various values of τ ≥ 0 and using several values of margin ν ∈ [0.01, 0.1]. The key point is that we can obtain good numerical results "even if" profiles are finally very sharp, using the proposed continuation method, consisting of numerically drive υ → 0. We have experimented using the conjugate gradient algorithm of the GSL (http://www.gnu.org/s/gsl) library, but have also checked that this is not a critical choice. Robustness has been checked for γ ± t = 1, φ t (s) = |s| and several for generalized metric also. Further numerical results are provided as supplementary material of this submission, while the code is available in the open-source EnaS library (http://enas.gforge.inria.fr).

Conclusion

The key point, here, is the non-learnability of even-based networks [START_REF] Šíma | On the nonlearnability of a single spiking neuron[END_REF], since it is proved that this problem is NP-complete, when considering the estimation of both weights in the general case, except for exact simulation [START_REF] Rostro-Gonzalez | Backengineering of spiking neural networks parameters[END_REF]. We show that we can "elude" this caveat and propose an alternate efficient estimation mechanism, inspired by alignment metrics used in spike train analysis [START_REF] Houghton | Measuring representational distances -the spike-train metrics approach[END_REF], thus providing a complement of other estimation approaches [START_REF] Schrauwen | Towards Applicable Spiking Neural Networks[END_REF], beyond usual convolution metric [START_REF] Houghton | Measuring representational distances -the spike-train metrics approach[END_REF][START_REF] Bruno Cessac | Overview of facts and issues about neural coding by spikes[END_REF]. At last, the proposed mollification is a series of convolution metric, but that converges towards the expected alignment metric.

A On modeling choices regarding deterministic discretetime models.

At the microscopic biological level, units are punctual models of neurons, events are action potentials, weights stand for synaptic strength, and the dynamical system is the neural network under consideration. At the mesoscopic biological level, units are cortical columns (see, e.g., [START_REF] Chemla | Biophysical cortical column model for optical signal analysis[END_REF] for a discussion on the concept), events are synchronization, rhythms, or sudden activity change, weights are related to the average connection strength, and the dynamical system is a cortical map (see e.g. [START_REF] Thierry Viéville | How do high-level specifications of the brain relate to variational approaches?[END_REF] for a discussion on the concept), including dynamical neural fields (see, e.g., [START_REF] Rougier | Dynamic neural field with local inhibition[END_REF]).

Considering discretized temporal mapping.

Considering that time is discretized is a twofold issue. On one hand, it corresponds to the fact that not all continuous time sequences correspond to event trains, since they are constrained by the network dynamics, yielding global time constraints such as the fact that inter-event intervals are bounded by a refractory period r and the fact that event times are defined up to some absolute precision δt (see [START_REF] Bruno Cessac | Overview of facts and issues about neural coding by spikes[END_REF] for a detailed discussion), this being true in both biological and electronical implementations. The maximal amount of information, for one unit, is thus bounded during a finite period [0, T ] as stated in [START_REF] Bruno Cessac | Overview of facts and issues about neural coding by spikes[END_REF]:

T r log 2 T δt bits.
In a biological context, the order of magnitude is of 1K bits/second for a neuron, in coherence with biological observations [START_REF] Rieke | Spikes, Exploring the Neural Code[END_REF]. On the other hand, at a pragmatic level, time discretized network models are rather easy to study theoretically [START_REF] Cessac | A discrete time neural network model with spiking neurons. rigorous results on the spontaneous dynamics[END_REF][START_REF] Cessac | On dynamics of integrate-and-fire neural networks with adaptive conductances[END_REF], trivial to simulate (contrary to continuous time models, see [START_REF] Brette | Simulation of networks of spiking neurons: a review of tools and strategies[END_REF] for discussion), and correspond without bias to what happens in a computer. We thus focus on discrete time and are going to briefly point out, from step to step, to which extent the present development can be applied to continuous time frameworks.

Considering a deterministic framework.

We also consider here a deterministic framework, i.e. work in a context where not the "average" input/output response, but an exact or approximate specific input/output response is targeted (i.e. in a context where "each spike may matter" [START_REF] Guyonneau | Neurons tune to the earliest spikes through stdp[END_REF][START_REF] Delorme | Network of integrate-and-fire neurons using rank order coding b: spike timing dependant plasticity and emergence of orientation selectivity[END_REF], which seems to be the case, e.g., in the biological visual system working in natural scenario [START_REF] Baudot | Nature is the code: high temporal precision and low noise in V1[END_REF][START_REF] Masland | The unsolved mystery of vision[END_REF][START_REF] Koch | Efficiency of information transmission by retinal ganglion cells[END_REF]). More precisely, from [START_REF] Bruno Cessac | Overview of facts and issues about neural coding by spikes[END_REF], we propose the following pragmatic view of the network result coding scheme (i.e. the "neural code" in a biological context [START_REF] Vanvreeswijk | What is the neural code? 23 Problems in System neuroscience[END_REF][START_REF] Rieke | Spikes, Exploring the Neural Code[END_REF]): two results correspond approximately to the same code if their distance with respect to a given metric or pseudo-metric 1 is small. For instance, considering a rank coding scheme (i.e. in a context where the relative temporal sorting of the events matter, but not their exact temporal values [START_REF] Delorme | Network of integrate-and-fire neurons using rank order coding b: spike timing dependant plasticity and emergence of orientation selectivity[END_REF]) the related pseudo-metric is discrete and easy to state: network output are equivalent if the ranks of event trains match, and not-equivalent otherwise. Contrary to this binary choice, our proposal is to introduce a richer structure: The proposed modeling view is not only to consider a weak notion of network coding where two codes can only be either equal or different, but a more general notion where two codes are similar up to some quantified distance. This seems to correspond to a more realistic view of, for instance, the still mysterious "neural code" (see [START_REF] Bruno Cessac | Overview of facts and issues about neural coding by spikes[END_REF] for a discussion), and at the estimation level allows variational optimization mechanisms to be used.

B Variational formulation of the alignment metric

Let us derive and explain the choice of (2). We proceed in three steps.

Coincidence metric mollification. We start with the simple case where τ = 0.

In that case we have no shift but only insertion/deletion. We count γ + t if we have to insert a spike (Z n [t] = 0 and Zn [t] = 1) and count γ - t if we have to delete a spike (Z n [t] = 1 and Zn [t] = 0). The total counts writes:

d 0 Z, V τ =0 def = nt γ ± t | Zn [t] -Z n [t]| = nt γ ± t H (1 -2 Zn [t]) (V n [t] -θ) , (3)
writing:

γ ± t def =    γ + t , Zn [t] = 1, δ[t] = 0 γ - t , Zn [t] = 0, δ[t] = 0 0 , δ[t] = 0. (4) This comes from the fact that | Zn [t] -Z n [t]| ∈ {0, 1} is equal to 1 if an only if Zn [t] = Z n [t]
, while the second form derives from the fact, obvious to derive,

1 A metric, corresponding to the intuitive notion of "distance" d(R, R ′ ) between two results R and R ′ , is a real function, positive (d(R, R ′ ) ≥ 0), symmetric (d(R, R ′ ) = d(R ′ , R)), definite (d(R, R ′ ) = 0 ⇔ R = R ′ ))), thus semi-definite (d(R, R) = 0), and subadditive (d(R, R ′ ) ≤ d(R, R ′′ ) + d(R ′′ , R ′ ), given any third result R ′′ ).
If such a function is positive, definite, subadditive, but not symmetric, it is called a divergence. If such a function is positive, symmetric, subadditive, but only semi-definite while not definite, it is called a pseudo-distance. Thanks to the subadditivity, a pseudo-distance induces a distance on the quotient space of the equivalence classes of elements at a zero distance to of each-others. On the reverse, an equivalence relation ≡, corresponds to a discrete pseudo-metric

(d(R, R ′ ) ∈ {0, 1}) such that R ≡ R ′ ⇔ d(R, R ′ ) = 0.
It is the coarser pseudo-metric compatible with the related equivalence classes, while finer distances quantify to which extents two results differ (see, e.g., [START_REF] Dieudonne | Treatise on Analysis[END_REF] for a text book reference).

that:

| Zn [t] -Z n [t]| = | Zn [t] -H (V n [t] -θ) | = H (1 -2 Zn [t]) (V n [t] -θ) .
If γ ± t = 1 the distance is the Hamming distance, counting the number of non-coincidences.

This distance is discontinuous, with a jump each time V n [t] crosses the threshold θ. We now replace the Heaviside function by a mollification, H υ in (3). The exact choice of H υ is detailed in section C. The corresponding mollified metric d υ ( Z, V), is now explicitly varying with V, not only the value of Z, allowing to tune the state value in order to it drive towards the correct side of the threshold.

The term

H υ (1 -2 Zn [t]) (V n [t] -θ) now increases the value of V n [t] if Zn [t] = 1 and decreases the value of V n [t] if Zn [t] =
0 and, thus drive the the state value towards the correct threshold side.

The interest of using the second form of ( 3), is that it has a symmetric effect for values below and above the threshold, which is not the case with the former form.

Variational form of the alignment metric. We now consider an alignment metric with its corresponding indexing δ.

As for the distance calculation or the divergence indexation, we can construct the variational form by double induction, on the n-th value such that Z n [t n ] = 1 and the n-th value such that Zn [ tn ] = 1. This is due to the fact that, in a minimal path, each event can be either deleted or shifted once to coincide with an event in the other event-train. Also, an event can be inserted only at a time that matches the occurrence of an event in the other event-train.

Furthermore, a minimal path cannot include an insertion of an event that is later deleted or shifted, or a deletion of an event that is later inserted or shifted, or a shift in both direction, since the cost of such path can be reduced by eliminating some steps. Individual events cannot intersect, since uncrossing them reduces the amount of shifting. This means that t + δ[t] is an increasing function, i.e.,:

∀t ′ < t, Zn [t ′ + δ[t ′ ]] = Zn [t + δ[t]] = 1 ⇒ t ′ + δ[t ′ ] < t + δ[t] ⇔ -1 < δ[t]-δ[t ′ ] t-t ′ . (5)
We write {t 1 , t 2 , • • • }, i.e., t k , k > 0 the k-th value such that Z[t k ] = 1, with a similar notation for tk, k > 0.

We start with dk ,0 = l<k γ - tk due to the fact that the distance between any event-train and the empty event-train corresponds to the cost of deleting all events, while δ(tk) = -0 in this case. Similarly, d 0,k = l<k γ + t k corresponds to inserting all events, with δ(tk) = +0 in this case.

In order to calculate dk +1,k+1 we may either perform a deletion from dk ,k+1 , an insertion from dk +1,n or a shift from dk ,n as made explicit in (1).

-If there is deletion, we obtain dk +1,k+1 = dk ,k+1 + γ - tk , while Z n [t k+1 ] = 0 and Zn [ tk] = 1. And we must have Z n [ tk] = 0 also, otherwise it would be costless not to perform the deletion. The cost to take into account is γ - tk since the event deletion occurs at time tk.

-If there is an insertion, in a similar way, dk +1,k+1 = dk +1,n + γ + t k , while Z n [t k ] = 1 and Zn [t k ] = 0.

-If there is shift, we obtain dk +1,k+1 = dk ,k +φ t k tk-tk τ

, with Z n [t k ] = 1 and Zn [ tk ] = 1, while we must have Zn [t k ] = 0 (otherwise a trivla better alignment occurs at t k ).

All together, we obtain the induction given in [START_REF] Rieke | Spikes, Exploring the Neural Code[END_REF]. Now if we want to derive a form similar to (3) we must recalculate, given a function δ, i.e., given the choices between insertion/deletion and shifts being fixed, the alignment function, now from an induction over the time t. To this aim, we write d {0,t{ = d Z, Z 0≤t ′ <t . For t = 0, obviously d {0,0{ = 0. Assuming we have calculated d {0,t{ at time t, we may have the following cases: . The key point is that we claim these are the only possible cases: -when (0, 0) nothing happens, -when (1, 0) the event Z n [t] can only be inserted or the source of a shift, -when (0, 1) the event Zn [t] can only be deleted or the target of a shift, -when (1, 1) it may either be a coincidence of both events, or the result of two shifts, other operations being not optimal as discussed before. Furthermore, since the shift cost must be counted only once and is taken at the time t, it must be counted at the source but not the target of the shift.

Z n [t], Zn [t] d {0,t+1{ -d {0,t{ (0, 0) 0 no event (1, 0) γ + t insertion of Zn [t] (1, 0) φ t shift source from to Z n [t] to Zn [t + δ[t]] (0, 1) γ - t deletion of Zn [t] (0, 1) 0 shift target at Zn [t] from Z n [t ′ ], t = t ′ + δ[t ′ ] (1, 1) 0 coincidence (1,
All this can be summarized writing:

d {0,t+1{ = d {0,t{ + γ ± t | Zn [t] -Z n [t]| + φ t
with γ ± t defined from (4) and it is obvious to verify that the formula is true for the seven possible cases.

As a consequence, for a event train Zn [t], t ∈ {0, T { and given an optimal indexing δ, we can write the corresponding alignment metric as d Z, Z = distance corresponds to the minimal value among all indexing, i.e. we can write 2 :

δ = argmin δ d {0,T { , d {0,T { = t γ ± t | Zn [t] -Z n [t]| + φ t .
Adding Kuhn-Tucker multipliers µ t allows to guaranty, using [START_REF] Bruno Cessac | Overview of facts and issues about neural coding by spikes[END_REF], that δ is an increasing function:

δ = argmin δ d {0,T { + t,Z[t]=1,t ′ =argmin t ′ >t Z[t ′ ]=1 µ t 1 + δ[t ′ ] -δ[t] t ′ -t .
The minimization is to be initialized with δ[t] = 0, i.e. considering only deletion/insertion and no shift.

A step further, in order to solve the ambiguity between indexing functions, from the last time to the previous time, we must choose a function γ ± t that decreases in the past, while its value must bounded by the higher shift cost, i.e. φ t < γ ± t , for the maximal acceptable value of δ[t].

Mollification of the alignment metric. Let us finally proposes a mollification derived from the mollification proposed for the coincidence metric and taking into account the shifts defined by the alignment metric. Here we consider that the δ indexation is fixed, and propose as criterion:

d υ Z, Z = nt υ d {0,T { , υ d {0,T { = γ ± t H υ (1 -2 Zn [t + δ[t]]) (V n [t] -θ) + φ t Z n [t] + (1 -Z n [t]) H υ -δ[t] τ H υ (θ -V n [t]) ν=0 .
Obviously, the term related to γ ± t enjoys the same properties as for the coincidence metric and we simply re-use what has been developed in this case.

The term related to φ t has a different behavior depending on Z n [t]]: -If Z n [t] = 1, i.e., when there are both a shift source and a shift target at time t, there is no adjustment on V n [t] since reducing one shift delay may increase the other one, with an unpredictable effect.

-If Z n [t] = 0, i.e., when there is only a shift source at this time t, then the value of V n [t] is increased/decreased depending on the sign of δ[t] in order to reduce the absolute value of the delay:

δ[t] > 0 t < t better shift t in the future V n [t] ↓ δ[t] < 0 t < t better shift t in the past V n [t]
↑ in accordance with the signs of variation in the criterion. Furthermore, the variation of V n [t] is bounded in order the adjustment not to generate an additional 2 Interesting enough is the fact we can apparently easily define a continuous form of this variational definition writing:

d (z, z) = min δ R t γ(t) ± h |zn(t + δ(t)) -zn(t)| + φt " δ(t) τ "i + µ(t)(1 + δ ′ (t))
with δ(t) = 0 as initial value, writing z(t) = P spurious event. This is the reason why we consider

H υ (θ -V n [t]) ν=0 in order to maintain θ > V n [t].
This last part if the specification is clearly a heuristic that has been adjusted when experimenting at the numerical level, in order to observe the proper behavior.

Conclusion.

We have made explicit the different elements allowing us to propose (2) as a reasonable criterion. Given δ, the weights are adjusted in order to optimize the adequacy between the desired and expected event trains. At the implementation level, a relaxation minimization scheme is proposed:

1. Given a network event dynamic Z calculates δ to estimate the alignment distance;

2. For this indexing δ optimize the weights to tune V n [t] in order to reduce the alignment distance;

3. As soon as the event train is changed, repeat step 1 in order to re-estimate δ, stopping when a local minimum is found.

C Mollification of the Heaviside function

Let us specify how to mollify the event generation function ρ() = H(u -θ). The Heaviside function is defined here with H(0) = 0. In other words, we focus on the fact an event is defined by a state value above a given threshold θ. The generalization to other semi-algebraic conditions is straightforward, since they always can be stated as a combination of Heaviside functions (see [START_REF] Benedetti | Real algebraic and semi-algebraic sets[END_REF] for a treatise on the subject).

In words, we precisely need to replace H(u) by a regular function that can influence the estimation, -either if the condition is incorrect, -or if the condition is correct, but close to be incorrect, i.e. at margin ν of the correctness boundary, while -we better require the function to have no influence if the condition is correct and beyond this boundary.

A suitable function that fits with this requirement is the H υ (u) non-linear regular profile described in Fig. 3, and is defined as:

H υ (u) def = H(u + √ υ ν) exp - υ u + √ υ ν , (6)
thus with

H υ (u) = 0, u ≤ - √ υ ν and H υ (+∞) = 1.
Here ν is a margin, i.e. it allows one to maintain the state at a non-infinitesimal distance to the threshold. It is easy to verify that this function is regular including in 0.

It is also important to notice that the function is convex on ] -∞, υ 2 -√ υ ν[, i.e., when H υ (u) < e -2 (drawn in magenta in Fig. 3).

Furthermore, lim υ→0 H υ (u) = H(u), assuming that

H(0) = 0. This conver- gence is uniform in R -[- √ υ ν, - √ υ ν + a[, a > 0, thus in [0, +∞[. This comes from the fact that |H(u) -H υ (u)| < (1 -e - √ υ
), thus uniformly bounded when υ → 0. This would not have been the case without a margin ν > 0. More precisely, the reader can verify, that this is the denominator with an additive term υ d ν, d < 1, say d = 1/2, which guaranties the convergence and the uniformity of this convergence in [0, +∞[. Around -√ υ ν, however, it is not possible to have convergence uniformity, which seems to be intrinsically related to the discontinuity of H. As a consquence, at the numerical level values must not remain in this neighborhood, to avoid "spurious jumps" when υ → 0.

A step further, |H -H υ | 2 is bounded for the L 2 norm (which is not the case for the

L 1 norm), since +∞ -∞ (H(u) -H υ (u)) 2 du υ=1,ν=0 = 2 log(2), so that |H -H υ | 2 is
finite, which is also an important fact for numerical algorithms.

Since H υ (u) = 0, u ≤ -√ υ ν, while still regular, the influence of obviously correct events is avoided, allowing the criterion to better tune more critical events. More precisely :

H ′ υ (u) = H υ (u) υ (u + √ υ ν) 2 = H u + √ υ ν υ (u + √ υ ν) 2 + O 1 u , while H υ (u) < H(u) for u > 0.
The derivation is not obvious since formally:

H ′ υ (u) = b(u + √ υ ν) + H υ (u) υ (u+ √ υ ν) 2
with b(u) def = D (u) exp -υ u , D() being the Dirac distribution. However, for all compact domain bounded function φ(), we obtain: b(u) φ(u) = exp -υ u φ(u) u=0 = 0 thus b() = 0 in the distribution sense.

A step further, there is a closed-form for the n-th derivative, i.e.

H (n) υ (u) = H υ (u) n-1 k=0 (-1) k µ k,n υ n-k (u+ √ υ ν) 2 n-k = H υ (u) υ n (u+ √ υ ν) 2 n + O 1 u 2 n-1
, for positive integers, iteratively defined as:

µ k,n =    1 k = 0 µ k,n-1 + (2 n -k + 1) µ k-1,n-1 0 < k < n 0 n ≤ k
, straightforward to derive from a piece of computer algebra.

Morally, this mollification corresponds to a convolution H υ = υ * H with respect to some suitable convolution kernel (without, up to our best knowledge, a closed form expression).

All together, the choice of the mollification kernel is not trivial. Better choices do exists, but they have to enjoy at least the properties made explicit here.

D Numerical verification

In order to estimate the variational method performance in terms of precision, we consider the simplest possible linear feed-forward network:

V n [t] = m W nm I m [t] + β n [t], Z n [t] = H(V n [t] -θ)
where the the input I m [t] is a known deterministic signal and the additive noise β n [t] is drawn from a zero-mean random Gaussian distribution of standard deviation β. Using this model, we perform a master-slave validation test. The weights of the master are drawn from a Gaussian distribution of standard deviation σ and the additive noise of magnitude β, is added to the master output before presented to the estimation. The estimation is done on a slave noiseless linear feed-forward network of the same number of output units N and input units M , considering C trials of T time steps, which initial weights are similarly randomly drawn. Here we use N = 4, M = 2, T = 1000, C = 2 with θ = 0.5. And typically consider σ = 1, β = 0.25 unless specified. Default alignment metric parameters are υ = 1, µ = 0.1 and τ = 0 unless their values are explored.

Typical results are illustrated in Fig. 4, where it is numerically shown that the estimation using the proposed alignment variational estimation has the same performance as the "ideal" least-square estimation on the state values for small data errors. This is not a trivial result since the alignment variational estimation mechanism can not observed the master state values but only the related events. Furthermore, we approximate here a non-regular discontinuous criterion by the mollification mechanism and there were no guaranty that it should work. Here it does for small and reasonable errors (up to 30% of noisy events in this particular case). Fig. 4 also shows that the dispersion for this numerical test is small, yielding a robust result.

A step further, we numerically verify in Fig. 5 that the mollification parameters values have only a marginal influence on the final estimation, in this case. The use of a margin yields better results, as hypothesized in the specification, while the mollification introduces a small additional bias reduced when the mollified approximation converges towards the true alignment metric, these effects being of the seconder order. We also verify in Fig. 5 that these results are stable for different values of time-shift cost.

E Discussion.

Let us propose a short discussion about general links between the present development and what is known regarding the computational power of artificial or biologically plausible neural networks with event-based (i.e., spiking) activity [START_REF] Rieke | Spikes, Exploring the Neural Code[END_REF][START_REF] Maass | On the relevance of time in neural computation and learning[END_REF][START_REF] Dayan | Theoretical Neuroscience : Computational and Mathematical Modeling of Neural Systems[END_REF][START_REF] Gerstner | Spiking Neuron Models[END_REF]. Fig. 4: Estimation error bounds : Alignment error obtained without any estimation (in gray), with a least-square estimation on the state values (in red), and using the proposed alignment variational estimation (in green), for different additive noise β in abscissa. The ordinate corresponds to the average alignment error form 0 for an exact estimation to 0.5 for a totally random result. The alignment mollification estimation is bounded by the non-estimation case from above an ideal least-square estimation where all state values (and not only the events) could be known from below. Error dispersion : Representation of the numerical test dispersion for the alignment metric estimation. Error-bars correspond to ± one standard-deviation. The ordinate corresponds to the average alignment error. Time-shift cost error dependency : Comparison between the alignment error obtained without any estimation (in gray) with the proposed alignment variational estimation (in green), using the same conventions as for Fig. 4, showing the stability of the estimation not only for an event-count metric.

The key problem of calculability: non-learnability.

It is known that recurrent neuron networks with frequency rates are universal approximators [START_REF] Schäfer | Recurrent neural networks are universal approximators[END_REF], as multilayer feed-forward networks are [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF]. This means that neuron networks are able to simulate dynamical systems, not only to approximate measurable functions on a compact domain, as originally stated (see, e.g., [START_REF] Schäfer | Recurrent neural networks are universal approximators[END_REF] for a detailed introduction on these notions).

Spiking neuron networks have been proved to be also universal approximators [START_REF] Maass | On the relevance of time in neural computation and learning[END_REF] and that, theoretically, spiking neurons can perform very powerful computations with precise event timings. Spiking neurons are at least as computationally powerful as the sigmoidal neurons traditionally used in artificial neuron networks [START_REF] Maass | Fast sigmoidal networks via spiking neurons[END_REF][START_REF] Maass | Networks of spiking neurons can emulate arbitrary hopfield nets in temporal coding[END_REF]. This result has been shown using a spike-response model (see [START_REF]Pulsed Neural Networks[END_REF] for a review) and considering piece-wise linear approximations of the membrane potential profiles. In this context, analog inputs and outputs are encoded by temporal latencies of event firings. It has been shown that any feed-forward (multi-layer) or recurrent analog neuronal network (e.g. Hopfield network) can be simulated arbitrarily closely by an insignificantly larger network of spiking neurons. The assertion holds even in the presence of noise [START_REF] Maass | Fast sigmoidal networks via spiking neurons[END_REF][START_REF] Maass | Networks of spiking neurons can emulate arbitrary hopfield nets in temporal coding[END_REF]. Such theoretical results highly motivate the use of spiking neuron networks for modelling and simulation purpose.

The key point, however, is the non-learnability of spiking neurons [START_REF] Šíma | On the nonlearnability of a single spiking neuron[END_REF], since it is proved that this problem is NP-complete, when considering the estimation of both weights and delays. Here we show that we can "elude" this caveat and propose an alternate efficient estimation, inspired by biological models.

We also have to notice, that the same restriction apply not only to simulation but, as far as this model is biologically plausible, also holds at the biological level. It is thus an issue to wonder if, in biological neuron networks, delays are really estimated during learning processes, or if a weaker form of weight adaptation, as developed here, is considered.

As far as this contribution is concerned, we consider a weak notion of biological plausibility: A simulation is biologically plausible if it verifies an explicit set of constraints observed in biology. More precisely, we have taken time constraints, shared by all dynamics, further called "general time constraints". The time constraints are based on biological temporal limits and appear to be very precious quantitative elements, both for estimating the coding capacity of a system and for improving simulations.

Considering such learnability constraints, i.e., how can artificial or biological systems by-pass such computational barrier ? By-passing the non-learnability barrier.

As pointed out previously, the weights estimation is proved to be NP-complete. This means that in order to "learn" the proper parameters we have to "try all possible combinations of delays". This is intuitively due to the fact that each delay has no "smooth" effect on the dynamics but may change the whole dynamics in an unpredictable way. This is the way proposed to elude this NP-complete problem by considering another estimation problem. Here we do not estimate one delay (for each synapse) but consider connection weights at several delays and then estimate a balancing of their relative contribution. This means that we consider a weak delay estimation problem.

The alternative approach is to estimate delayed weights, i.e. a quantitative weight value W ′′ ojd or W ′ ijd at each delay d ∈ {1, D}. Obviously, the case where there is a weight W ij with a corresponding delay d ij ∈ {0, D} is a particular case of considering several delayed weights, since we can write W ijd = W ij δ(d -d ij ), δ() being the Kronecker symbol in this case. In other words, with our weaker model, we are still able to estimate a neuron network with adjustable synaptic delays.

We thus do not restrain the neuron network model by changing the problem, but enlarge it. In fact, the present estimation provides a smooth approximation of the previous NP-complete problem.

The key idea of finding an approximate solution to the previous NP-complete problem, is instantiated at the implementation level, using mollified metrics, and the interest of such approach has been stated considering numerical examples.

The main limit of the present approach: polychronization.

A spiking network can polychronize, i.e., exhibit reproducible time-locked but not synchronous firing patterns within 1 millisecond precision. Polychronization can be viewd as a generalization of the notions of synchronization and synfire chains. Due to the interplay between the delays and a form synaptic plasticity, the spiking neurons spontaneously self-organize into groups and generate patterns of stereotypical polychronous activity.

In [START_REF] Izhikevich | Polychronization: Computation with spikes[END_REF], it has been shown that the number of co-existing polychronous groups far exceeds the number of neurons in the network, resulting in an unprecedented memory capacity of the system. The author speculates on the significance of polychrony to the theory of neuronal group selection and cognitive neural computations.

In [START_REF] Paugam-Moisy | Delay learning and polychronization for reservoir computing[END_REF], the network processing and the resulting performance is explained by the concept of polychronization, The model emphasizes that polychronization can be used as a tool for exploiting the computational power of synaptic delays and for monitoring the topology and activity of a spiking neuron network.

Taking such complex aspects of the neural code into account cannot be performed by any available metrics. New metrics, taking long term interactions into account have to be developed and this is a challenging issue. This is not the case here and a challenging perspective of the present development.

Fig. 2 :

 2 Fig. 2: Solving ambiguous equal distance alignments, in the case where a shift cost equals the insertion/deletion cost.

Fig. 3 :

 3 Fig. 3: Defining the mollification of the Heaviside function H().

  t n k ), D being the Dirac distribution, and also define the related mollification, as proposed here in the discrete case.

Fig. 5 :

 5 Fig. 5: Mollification parameters error dependency : Marginal influence on the final estimation of the mollification parameters υ in abscissa, for three margin parameter values µ = {0, 0.1, 0.2} in {green, magenta, cyan} respectively.The ordinate corresponds to the average alignment error. Time-shift cost error dependency : Comparison between the alignment error obtained without any estimation (in gray) with the proposed alignment variational estimation (in green), using the same conventions as for Fig.4, showing the stability of the estimation not only for an event-count metric.

n d {0,T { , using the previous formula. On the reverse, consider d {0,T { as a function of the indexing δ, the alignment

* Partially supported by the ANR KEOpS project and CORTINA associated team.