Noisy Laplace deconvolution with error in the operator

We adress the problem of Laplace deconvolution with random noise in a regression framework. The time set is not considered to be fixed, but grows with the number of observation points. Moreover, the convolution kernel is unknown, and accessible only through experimental noise. We make use of a recent procedure of estimation which couples a Galerkin projection of the operator on Laguerre functions, with a threshold performed both on the operator and the observed signal. We establish the minimax optimality of our procedure under the squared loss error, when the smoothness of the signal is measured in a Laguerre-Sobolev sense and the kernel satisfies fair blurring assumptions. It is important to stress that the resulting process is adaptive with regard both to the target function's smoothness and to the kernel's blurring properties. We end this paper with a numerical study emphazising the good practical performances of the procedure on concrete examples.

Introduction
Laplace deconvolution is motivated by a wide set of practical applications, ranging from population dynamics or physics to computational tomography or fluorescence spectroscopy (Linz [22,Chap. 2], Ameloot et al. [4], Comte et al. [9]). In the corresponding setting we observe q, the result of the action of a kernel g on the function of interest f , according to the following equation Equation (1) is also refered to as Volterra integral equation. One of its main features is its causal property, since qptq is affected only by the values of f and g at times anterior to t. Of course, only finite samples of qptq are accessible in practice. Moreover, the presence of additional noise justifies the empirical modelization of (1) by the classical regression model, inspired by Abramovich et al. [3] ypt i q " ż t i 0 gpt i´τ qf pτ qdτ`ση i , i " 1, ..., n where 0 ď t 1 ď ... ď t n ď T n are the points of observation, pη i q i"1,...,n are independent standard gaussian variables, and σ is a fixed factor accounting for the precision of the observations. T n is supposed to grow with the number of observations n.
As pointed out in Abramovich et al. [3] and Comte et al. [9], in spite of its apparent similarity with the Fourier deconvolution problem, the theoritical features of equation (1), as well as the practical problems raised during its resolution are deeply different. More precisely, setting artificially gptq " f ptq " 0 for t ă 0 amounts to solving the classical Fourier deconvolution problem ypt i q " ż Tn 0 gpt i´τ qf pτ qdτ`ση i , i " 1, ..., n A first notable objection is that the framework of classical Fourier deconvolution assumes periodicity of the function f and the kernel g on r0, T s, a meaningless notion when applied to a varying time set r0, T n s. Even more problematic is the fact that this modelization totally ignores the causal feature of Laplace convolution, creating unwanted interferences between different time sets. To finish, the manipulation consisting in artificially expanding q and g for t ă 0 creates artifacts on the estimated function at times t ă 0 as well.
Another approach is to treat equation (2) as a general ill-posed problem and apply a Tikhonov regularization (Golubev [15]). However the direct implementation of this method also destroys the causal nature of equation (1), and tends to oversmooth the solution (Cinzori and Lamm [7]). Subsequent adaptations which remedy these shortcomings are present in Lamm [20] and Cinzori and Lamm [7]. However in these works the time set is considered to be fixed. A more suitable theoritical tool in solving (1) is the use of Laplace transform, which allows to derive a closed form of the solution. However, its direct implementation is compromised by numerical problems, since the generic expression of the inverse Laplace transform is not easily computable in general. This motivates the widespread use of inversion tables, unfortunately irrelevant when the image function is not known exactly but approximated via a numerical scheme. In this paper, following Comte et al. [9], we will exploit the properties of Laguerre functions, which can be used either to compute the inverse Laplace transform (Abate et al. [1], Lien et al. [21]), or to solve directly equation (1) (Keilson and Nunn [19]). More precisely, a Galerkin method applied to (2) shows that, even if their role is not entirely symmetric to the role played by harmonics in the framework of Fourier deconvolution, they allow a sparse analysis of equation (1). All the previous mentionned works only concerned the case of a deterministic noise at best. The presence of random noise requires an additionnal treatment, and calls for specific statistical tools. In the setting of random noise, Dey et al. [11] considered a kernel of the form e´a t and used a regularized inversion of the inverse Laplace transform. More recently, Abramovich et al. [3] conceived an optimal procedure in the minimax sense on Hölder spaces H s pR`q. This procedure used an exact expression of the solution involving the derivatives of q, which were then estimated via Lepskii's method. However a shortcoming of the procedure is its strong dependence on the kernel g, in the sense that a small error in g can translate into a wide difference in the result. In other words there seems to be a trade off between the closed form of the solution, and the unstability with regard to the kernel. Moreover, the fact that g is seldom observed directly in practice, but is usually subject to experimental noise should prompt us to privilege stability over exactitude.
In that spirit, Comte et al. [9] took advantage of the algebraic properties of Laguerre functions in the context of (2). With an adequate penalty term, they proposed an estimator which mimicks the oracle risk to within logarithmic terms. This modelization has the non negligible advantage of practical simplicity and efficiency, since solving equation (1) amounts to the inversion of a lower triangular Toeplitz matrix. Even if this latter procedure proves to be more stable with regard to g experimentally, no systematic study has been conducted on the subject yet. In this paper we attempt to fill in this gap: we suppose that the observation of g is contaminated by a gaussian white noise, and show how Laguerre functions allow to handle this issue. We place ourselves under the minimax point of view and suppose that f belongs to a Laguerre-Sobolev space and that g satisfies standard blurring assumptions. We apply recent techniques for the treatment of noisy operators in the context of inverse problems (Hoffmann and Reiß [17],Delattre et al. [10]), which consist in a preliminary processing of the operator K coupled with a classical thresholding procedure applied to y.

Laguerre functions
Suppose that the target function f and the kernel K both lie in L 2 pR`q. Define the Laguerre polynomials (see Gradshteyn and Ryzhik [16]) and, following Comte et al. [9], the ensuing Laguerre functions, depending on the parameter a ą 0, The parameter a is a tuning parameter used to fit experimental curves. The Laguerre functions constitute a Hilbert basis of L 2 pR`q. Any function f P L 2 pR`q satisfies The following proposition illustrates the conveniency of Laguerre functions in the framework of equation (1). @a ą 0, @t ě 0, From now on, except if explicitly mentionned, we will suppose a " 1 2 .

Galerkin method
Proposition 2.1 prompted Comte et al. [9] to apply a Galerkin scheme to equation (1). Galerkin schemes rely on the choice of a set of functions which discretize the inverse problem at stake in a convenient way. They were beneficially applied in the context of inverse problems (Cohen et al. [8]), and blind deconvolution (Efromovich and Koltchinskii [14], Hoffmann and Reiß [17] and Delattre et al. [10]). To this end we will remind briefly the underlying methodology of a Galerkin scheme and show how it conveniently applies to equation (1). Let f P L 2 pR`q and K an operator of L 2 pR`q, and suppose we want to recover f from the observation q " Kf . Note V ℓ the finite dimensional space spanned by the orthogonal set of Laguerre functions tϕ k u kďℓ . The Galerkin approximation f ℓ of f on V ℓ is the solution of the equation We shall note K ℓ the Galerkin matrix pK ℓ q i,j " xKϕ j , ϕ i y, i, j ď ℓ. Note hence K the operator of L 2 pR`q mapping f onto t Þ Ñ ş t 0 f pt´τ qgpτ qdτ . We can reformulate (1) as Moreover, Proposition 2.1 implies: Proposition 2.2 (Comte et al. [9], Lemma 1). The Galerkin matrix K ℓ is lower triangular, Toeplitz. More precisely, note 9 g the function with Laguerre coefficientš In the sequel, for any function f P L 2 pR`q, we will note T pf q the infinite Toeplitz matrix such that T pf q i,1 " f i`1 for all i ě 0, and T ℓ pf q the extracted matrix defined by T ℓ pf q i,1 " T pf q i,1 , i ď ℓ`1. In particular, The resolution of the linear system (8) now shows great practical conveniency, provided that K ℓ is invertible. This is equivalent toǧ 0 ‰ 0, an assumption we will make in the sequel.

Application to the regression model with irregular design
It remains to incorporate two supplementary features of equation (2) in the inversion of (9). First, the presence of the random noise η and secondly, the possible irregularity of the design points. This construction is due to Comte et al. [9]. Due to the fact that the observation points t i are imposed by the problem, the estimation of the Laguerre coefficientsq ℓ of the function q suffers from two potential drawbacks. First, the infinite support of the Laguerre polynomials as well as the function q which should not be too problematic, provided that T n is large enough and that the functions decrease sufficiently to infinity. More problematic is the fact that the observation points t i are sometimes subject to experimental constraints, which affect their repartition on R`. The consistency of the estimation ofq ℓ is hereby deteriorated. We will hence suppose that the following conditions are fulfilled: ‚ There exists an integer n 0 such that n Tn ą σ for all n ě n 0 . ‚ lim nÑ8 T n " 8, and lim nÑ8 T n n " 0 To take into account the irregularity of the design, we follow Comte et al. [9] and define P n : r0; T n s Ñ r0; T n s a regular non decreasing function such that P n p0q " 0, P n pT n q " T n , P n pt i q " i n T n for i ď n Note Φ ℓ the pℓ`1qˆn matrix with entries pΦ ℓ q k,i " ϕ k pt i q. For any function h P L 2 pRq, we have where P ℓ is the orthogonal projector onto V ℓ . We deduce that where η ℓ " N p0, I n q. Let us take a closer look to the matrix p t Φ ℓ Φ ℓ˘. Its general term is ż Tn 0 ϕ k pτ qϕ ℓ pτ qP 1 pτ qdτ for n, T n large enough. If the points t i are equispaced, taking P pτ q " τ in (10) entails that T n n´1p t Φ ℓ Φ ℓ˘i s close to the identity provided that T n is large enough. As in Comte et al. [9], we hence reformulate (11) as the sequential model , Ω ℓ q and Ω ℓ " nT´1 n p t Φ ℓ Φ ℓ˘´1 . In general, Ω ℓ somehow quantifies the distance to the uniform design case. To ensure that the design is not too ill conditionned, we will suppose that the following assumption is fulfilled.
Assumption 2.3. Let L P N. There exists C ě 0, such that for all ℓ ď L, for all λ P SppΩ ℓ q, λ ď C This assumption is dependent on the integer L, which plays the role of a maximal resolution level, and will be adapted to the case of interest later. The inversion of (11) now requires controls of the variable pK ℓ q´1ξ ℓ . Under suitable properties of f and g, we shall be able to apply a classical inverse/thresholding procedure, and derive rates of convergence over specific regularity spaces. These properties are the subject of Part 3.

Error in the operator
We already mentionned the fact that the resolution of (1) is usually unstable with respect to g (Abramovich et al. [3]). Furthermore, in practice, inference on the kernel g is possible only through experimental noise, and requires a preliminary step of estimation giving way to imprecision. This additionnal error might significantly contaminate the result of any procedure of estimation if not properly treated. Let us see how Laguerre functions ϕ ℓ allow to handle this issue: in section 2.2, we established that the discretization of (13) with Laguerre functions involved a Toeplitz matrix with entries constituted of the Laguerre coefficients of 9 g. We can thus consider 9 g as the finite impulse response of the operator K when applied to the system pϕ ℓ q ℓě0 . To take into account the imprecision in the observations of 9 g, we adopt the framework of blind deconvolution and suppose that 9 g is not known exactly, but that we have acces the noisy version 9 g δ " 9 g`δb (12) where b is a gaussian white noise on L 2 pR`q. The generic problem of blind deconvolution is motivated by numerous scientific fields, including for example electronic microscopy or astrophysics, where the corresponding kernel is seldom known nor directly observed. It was adequatly discussed in Efromovich and Koltchinskii [14] and Hoffmann and Reiß [17]. Taking into account the observations (12), the projection 9 g ℓ is changed to 9 g ℓ δ " 9 where B ℓ " T ℓ pbq is a random Toeplitz matrix. In the sequel, for the sake of clarity, we note ε " σ b Tn n . Remark 2.4. We could as well suppose that we observe g δ " g`δb, yet it is more convenient to work with 9 g (the entries of the noisy Toeplitz matrix B are directly i.i.d standard gaussian variables). In the former case, the rest of the paper however adapts with no change in the algorithms, since inequality (27) is satisfied as well. A modification of the proof of Theorem 4.6 should also provide the lower bound for the second procedure.

Sobolev spaces associated to Laguerre functions
We proceed to the description of regularity spaces associated with the resolution of (13). The following material is classical, we refer to Bongioanni and Torrea [5] or Rathnakumar [26] for example. Since f Þ Ñ ? 2af p2a.q is an isometry of L 2 pR`q, the structures defined for different values of a are equivalent. Hence we shall only concentrate on the mainstream case where a " 1{2. Define the operator L on L 2 pR`, dxq by The functions ϕ ℓ are the eigenfunctions of L associated with eigenvalues pℓ`1 2 q. We hence define the Sobolev space W s associated with Laguerre functions as For a function f P L 2 pR`, dxq, we have the straightforward equivalence xf, ϕ ℓ yˇˇ2 ă 8 and the associated norm For M ě 0, we shall note W s pMq the Sobolev ball of radius M. Finally, we remind that, as }ϕ ℓ } 8 ď 1 for all ℓ ě 0, we have s ą 1{2 ñ W s Ă C 0 pR`q. From now on, we will hence suppose that there exists s ą 1{2 such that f P W s .

Banded Toeplitz matrices
Before entering into details about the kernel features, we introduce basic material on Toeplitz matrices. Most of it is inspired by Böttcher and Grudsky [6] and Comte et al. [9]. Let a " pa ℓ q P ℓ 1 pZq be a sequence of real numbers. We remind from section 2.2 that we note T paq the infinite Toeplitz matrix defined by The Toeplitz matrices T paq and T ℓ paq are naturally linked to the two respective Laurent series apzq " We will indifferently refer to the vector a or the corresponding Laurent serie. The spectral norm of T paq is related to the behaviour of apzq, as illustrated in the following proposition.
Proposition 3.1. Let a P ℓ 1 pZq. Let C stand for the complex unit circle. We have a ℓ z ℓˇ. A simple corollary is the following inequality In particular, Proposition 3.1 applies to the case of truncated Toeplitz matrices T ℓ paq. Moreover, if a has no zero on the complex unit circle, we have lim sup ℓÑ8 }T ℓ paq} op ă 8 and lim ℓÑ8 }T ℓ paq} op " }T paq} op (15) Now suppose that a and a 1 both generate lower triangular Toeplitz matrices (i.e. a k " a 1 k " 0 if k ă 0). Then the following equalities hold for all ℓ ě 0: In other words, the matrix multiplication (resp. inversion) is equivalent to a power serie multiplication (resp. inversion).

Degree of ill posedness
We now need to precise the properties of K as a blurring operator of L 2 pR`q. Usually the operator K is not compact, and the problem (1) is ill-posed. This results in practical unstabilities when trying to invert equation (11) from discrete observations. The quantification of the ill-posedness of the problem is specified by the introduction of a constant, called degree of ill-posedness (DIP) of the problem (see Nussbaum and Pereverzev [25], Mathe and Pereverzev [23] for a generic review). We adapt this concept to our framework, and make the following assumption.
Assumption 3.2 (Degree of ill-posedness of g). There exists ν ě 0, Q ě 0 such that, for all ℓ ě 0, }pK ℓ q´1} op ď Qpℓ _ 1q ν ν is called degree of ill-posedness of g (or equivalently of K). We note K ν pQq the set of functions which satisfy this assumption.
We shall see examples of kernels satisfying this assumption further. For the moment, we concentrate on the treatment of observations (13) in the context we just described.

Algorithms and rates of convergence
The main challenge which remains to be treated now is to articulate the two critical steps of inversion and regularization, via adapted procedures. For example, let us give a brief overview of the methodology in Comte et al. [9]: Let ℓ P N, and let Λ be the following contrast function, defined on R ℓ by Λ : t Þ Ñ }t} 2´2 xt, pK ℓ q´1y ℓ y Note }.} op the spectral norm and }.} HS the Hilbert Schmidt norm. A model selection is performed on the maximal level L, by introducing the following penalizing factor (B ą 0 is an arbitrary constant): where Q ℓ " pK ℓ q´1Ω ℓ t pK ℓ q´1 and a Q ℓ is a lower triangular matrix satisfying a Q ℓ t a Q ℓ " Q ℓ . The maximal levelL is hence chosen as where ℓpnq is a large enough resolution level, possibly depending on n, and the ensuing estimator of f is pKLq´1yL We follow here a different path: we suppose that the target function belongs to a Sobolev-Laguerre space, and perform thresholding techniques in a minimax framework. Furthermore, our results are asymptotic with regard to ε, δ. Would g be known, the estimation of f from observations (13) amounts to solving a standard inverse problem with signal noise. To this end, a prolific litterature is at disposal (a selected list is Donoho [12], Abramovich and Silverman [2], Cohen et al. [8]). In order to take into account the presence of noise in the operator, we shall hence apply a preliminar regularizing thresholding procedure to the noisy operator K δ in order to ensure the stability of the further inversion step. To that end, define the maximal level as with λ a positive constant. Define also the two thresholding levels We call this procedure Algorithm I. The preliminary threshold performed on pK ℓ δ q´1 ensures its proximity with pK ℓ q´1 with high probability (see Lemma 6.2). We now study the squared loss performance of the procedure.
The rates in Theorem 3.3 reveal two components, accounting respectively for the imprecision in the observation of the operator and the signal. The latter is fairly classic in non parametric statistics (Nussbaum and Pereverzev [25], Johnstone et al. [18]) where it is also optimal, while the former is standard (and optimal too) in blind deconvolution on Hilbert spaces (Efromovich and Koltchinskii [14], Hoffmann and Reiß [17]). Thus, we do not study the optimality of these rates in this paper, but rather concentrate on a more specific framework related to the problem of interest.

Adaptation to the standard framework of Laplace deconvolution
We now discuss the adapation of our algorithm in the mainstream framework of Laplace deconvolution, as exposed in Abramovich et al. [3] or Comte et al. [9]. As we shall see, this more restrictive framework allows to treat observations (13) more efficiently. To this end, we first define a more restrictive version of the degree of ill-posedness. gq, ϕ k y, so that p1{ 9 gqpzq " ÿ kě0 γ k z k . There exists ν ą 0, there exists Q 2 , Q 1 ą 0, such that for all ℓ ě 0, For Q " pQ 1 , Q 2 q, we note G ν pQq the set of functions g P L 2 pR`q such that Assumption 4.1 holds. Note that the validity of this assumption automatically entails Q 1 ď`1`2 2ν 2ν˘Q 2 . Note also that the left term in (21) is the Hilbert-Schmidt norm of pK ℓ q´1. Thus, Assumption 4.1 is more restrictive that Assumption 3.2. However, it is satisfied by a natural class of functions g: Proposition 4.2 (Comte et al. [9], Lemma 3/ Lemma 5). Suppose that there exists C, ν ą 1{2, µ P C and wpzq " ś N i"1 pz´µ i q, |µ i | ą 1 a polynomial function with no pole inside of the complex unit disc, such that Then Assumption 4.1 is satisfied. Furthermore, if w " 1 and ν ě 0, then |γ ℓ | " ℓ ν´1 Γpνq . For completeness, we give a proof of Proposition 4.2 in section 6. We now turn to the standard framework of Laplace deconvolution, as exposed in Abramovich et al. [3] and Comte et al. [9]. To this end, we define the following assumptions concerning the kernel g. (A1) There exists an integer r ě 1 such that (A2) g P L 1`r 0,`8q˘is r times differentiable and g prq P L 1`r 0,`8q˘.
(A3) The Laplace transform of g has no zeros with non negative real parts except for the zeros of the form 8`ib.
The consequences of these assumptions are well formulated in the terms of the preceding framework: Hence, Assumption 4.1 is verified with ν " r and Algorithm I applies. However, Assumption 4.1 provides additional information on the behaviour of p1{ 9 gq. We adapt Algorithm I to this new framework, by operating the following changes: ‚ Set the maximal level to where }A} HS " a Trp t AAq is the Hilbert-Schmidt norm. We call the modified procedure Algorithm II and note r f II the corresponding estimated function. A notable gain of this new algorithm is its independence with regard to the parameter ν. Indeed, Assumption 4.1 allows us to use }pK ℓ δ q´1} HS in (19) as a substitute of ℓ ν , and to overesimate the 'true' maximal level L I . Its performances are exposed in the next theorem: Thus, in addition to the adaptivity over the parameter ν, the strengthening of Assumption 3.2 via (20) and (21) allows to improve on the rates of Theorem 3.3 with regard both to the operator and signal noise. Our next result shows that the rate achieved in Theorem 4.5 is indeed optimal, up to logarithmic terms. The lower bound will not decrease for increasing noise levels δ and ε, whence it suffices to provide separately the cases δ " 0 and ε " 0. Theorem 4.6. Let s ą 1{2, let M ě 0 ν ą 1{2 and Q 2 ě c ν Q 1 ą 0. Here c ν is a constant depending only on ν which will will not seek to precise. We have where the infimum is taken among all estimatorsf of f based on observations (13).
Combining Theorem 3.3 together with Theorem 4.6, we conclude that our algorithm is minimax over W s pMq to within logarithmic terms in ε and δ, uniformly with regard to the blurring kernel g P G ν pQq.

Practical performances
In this section we study the practical performances of the two procedures developped above. Note that three potential sources of errors may contaminate the quality of the observations in (13) : the signal precision σ b Tn n , the operator precision δ and the design quality }Ω ℓ } op . We shall hence emphasize their influence in the estimation of f , as well as their respective interactions. Our first aim is to study the interaction between the effect of signal and kernel noise in the two procedures of reconstruction. To this end, we will isolate them from the effect of the design, and suppose that the latter is ideally conditionned by setting Ω ℓ " I ℓ . Let us start by a few precisions concerning the tuning parameters of Algorithm I and II. The setting up of these procedures requires the preliminary definition of λ, κ, τ sig and τ op . Tuning parameters: for the definition of the maximal level of resolution, we set λ " 1 for both algorithms. The concrete choice of adequate thresholding constants κ and τ is a complex issue. Our practical choices will be based on the following remark, inspired by Donoho and Johnstone [13]: in the case of direct estimation on real line, the universal threshold which is both efficient and simple to implement, takes the form 2 a | log ε|. A consistent interpretation is to consider that this threshold should kill any pure noise signal. We will adapt this reasoning to the case of interest. Choice of κ : we use as a benchmark the case where g " 0. Given δ large enough, we define κ as the smallest value κ δ such that , for all ℓ ď 10, 1 t}pK ℓ δ q´1}opăO´1 δ,l u " 0. The results are reported in Table 1 and give κ " 0.3. Choice of τ sig and τ op : It is clear that the role of τ sig and τ op is to control the influence of the signal (resp. the operator) error. To choose τ sig (resp. τ op ), we therefore set ε sig ą δ sig ą 0 (resp. δ op ą ε op ą 0) large enough. We resort to the case f " 0 as a benchmark: we have xf , ϕ ℓ y " 0 for ℓ ě 1, consequently the observations xg ε sig , ϕ ℓ y, ℓ ě 0 are pure noise. We hence simulate K δ sig and, integrating the precedently computed value of κ, apply the procedure for increasing values of τ sig (resp. τ op q until all the computed coefficients xf i , ϕ ℓ y (i "I,II) are killed for ℓ ď 10. The results are reported in Table 2.
κ 0.1 0.2 0.3 N 3 1 0 Table 1. Choosing of κ. N is the average number, computed on a basis of 10 realizations, of levels ℓ ď 10 such that }pK ℓ δ q´1} op ă O´1 δ,l pκq. We have δ " 10´2.  Table 2. Choosing of τ . For pδ sig , ε sig q " pε op , δ op q " p10´2, 10´1q and each value of τ , we computed 10 times the described procedure and reported N i the average number of remaining Laguerre coefficients for Algorithm i. We now apply the two procedures to the case where f 1 ptq " pt 2´t q expp´tq and g " ϕ 0 (a graphical representation of these two functions is presented in Figure 1).
We have`1 { 9 g˘pzq " p1´zq´1 " ÿ ℓě0 z ℓ hence Assumptions 3.2 and 4.1 are both satisfied taking ν " 1. For several values of ε and δ, we report the corresponding squared loss, computed on a basis of 500 realisations with the use of Parseval's identity, in Table 3. The corresponding results are presented in Figure 2 for one particular realization of ξ, b. The results indicate that the transition on the two types of errors occur when δ is higher than ε, translating a prevailing effect of the signal noise ε over the operator error δ in practice. As Theorems 3.3 and 4.5 suggest, the second Algorithm overperforms the first in (almost) every case. Discussion on the design irregularity: to control the squared risk of the two  Table 3. Normalized mean squared error of the two procedures applied to the functions f 1 and g. The computations were performed using a monte carlo method on 500 realizations. procedures, one needs condition 2.3 to be fulfilled. If not, the eigenvalues of the matrix Ω L become potentially too large, and observations (11) are not conveniently treatable. In this case, it is preferable to lower the maximal level down to a point where }Ω L } op remains under control. To this end, we change the maximal level of the two respective procedures to where α is an arbitrary thresholding constant, set to 1.5 in the sequel. We now fix σ " δ " 10´2 and chose the design points t i as t i " 100i{n for n " 200, 250, 750 and 1000. Taking the same kernel g " ϕ 0 , and setting f 2 ptq " pt 1{2´t q expp´tq, we compare the performances of the new choice N i to the previous one L i , by computing the respective mean squared losses on a basis of 500 observations and report the result in  Table 4. Normalized mean squared error of the two procedures when the design is constituted of 200 equispaced points on the interval r0; 100s. We compare the performances of the two maximal resolution levels L i and N i for the parameters σ " δ " 10´2, g " ϕ 0 and f 2 ptq " pt 1{2´t q expp´tq. since L i is usually already smaller than N i . However, the gain is notable for Algorithm II when n ď 250. To illustrate this point, we plot in Figure 3 the corresponding results when n " 200.

Back to the regression model
We now turn back to the original model (2) to apply the two procedures. It is well known that this model is asymptotically equivalent to (11), in the sense that a fine  enough design will provide an estimation of the Laguerre coefficients with a negligible error when n Ñ 8. We work with f 3 pzq " p1´zq 1{2 , g " ϕ 0 , δ " 10´2, and suppose that the design is constituted of the points t i " i ÿ j"1 pstep`|X j |q where pX j q jďn is an i.i.d sequence of N p0, 10´2q variables. We observe the noisy values ypt i q " qpt i q`ση i where qpzq " p1´zq 3{2 , and compute the Laguerre coefficientsq ℓ via the approximatioň We apply the two procedures and present the results on Figure 4.

Proofs
In the sequel, for the sake of clarity, we suppose that τ sig " τ op ∆ " τ .

Proof of Proposition 4.2
Proof. We can restrict ourselves to the case where µ " 1. Proposition 16 applied to equality (22) entails As a consequence, and }pK ℓ q´1} HS ě C´1}T ℓ`w˘}´1 op }T ℓ`p 1´zq´ν˘} HS Since w is assumed to have no zeros on C, we deduce from Proposition 3.1 that both }w´1} circ and }w} circ are finite, and from (15) that It remains to treat the binomial serie p1´zq´ν. This serie can be expanded as Γpℓ`1qΓp´ν´ℓ`1q is the generalized binomial coefficient. Furthermore, we haveˆ´ν ℓ˙" ℓÑ8 p´1q ℓ Γpνqℓ´ν`1 (26) which is a direct consequence of Euler's definition of the Gamma function Γpzq " lim kÑ8 k!k z Π k i"0 pz`iq . Since ν ą 1{2, the serie ÿ kˆ´ν k˙2 is hence divergent, and there exists r Q 2 , r Q 1 ą 0 such that, for all ℓ ě 0, The proof is complete thanks to (24) and (25).

Preliminary lemmas
We begin with the following lemmas. Lemma 6.1 is a concentration inequality on the variable }B ℓ } op , which results from a concentration inequality on subgaussian processes. Lemma 6.2 states that }pK ℓ δ q´1} op behaves as }pK ℓ q´1} op on a set with large probability. Finally, Lemma 6.3 establishes deviations bounds on the variables ζ ℓ´f ℓ which will be useful throughout the proofs of Theorem 3.3 and Theorem 4.5.
Lemma 6.1. There exists β 0 , c 0 independent from ℓ ě 0, such that, for all ℓ ě 0, for all t ě β 0 , This readily entails the following moments control, available for all ℓ ě 0, p ě 1 Proof. The proof is a slight modification of Meckes [24,Theorem 1], to which we refer for a complete study. Lemma 6.1 is trivially satisfied if ℓ " 0, 1, hence we will suppose that ℓ ě 2. From Proposition 3.1, we derive that We claim the two following facts: ‚ Let a 0 , ..., a ℓ P R. There exists c ě 0 such that ofor all t ą 0, The first point is readily verified since pb k q kďℓ is a standard Gaussian vector, while the second point directly results from the bounďˇe 2iπkx´e2iπkyˇď 2^2πk|x´y| for all x, y P r0, 1s, k ě 0 (see Meckes [24] for the rest of the proof). The deviation bound is now a consequence of Talagrand [27,Lemma 5.3]. Indeed, for all x P r0, 1s, which ends the proof.

Proof of theorem 3.3
Lemma 6.3. Under Assumption 3.2, we have, for all ℓ ě 0, Proof. In order to prove Inequalities (33) and (34), it suffices to study the tails of the random variablesˇˇxpK ℓ δ q´11 A ℓ 1 B ℓ`´δ B ℓ f ℓ`ε ξ L I˘, ϕ ℓ yˇˇ. For convenience we will only treat the case where ℓ ě 2, otherwise the result follows by identical arguments. To this end, we study each term apart. On A ℓ X B ℓ , Lemma 6.2 and Assumption 3.2 entail Thus, combining Assumption 2.3 with the latter inequality, a brief conditionning argument readily yields et us study the second term. On A ℓ X B ℓ , we have Let's now bound separately r 1 and r 2 . We first apply equality (16) to get where pb ℓ q 1 k " pb ℓ q ℓ´k . The result is a centred gaussian variable with variance which hence satisfies et us study the term r 2 . Since the maximal level L verifies L ď λpδ| log δ|q´1 ν`1 , we have, for all ℓ ď L, δℓ ν`1 log ℓ À 1. We deduce that À expp´t`δℓ ν˘´1 q1 ttąβ 0 δℓ ν u`1ttďβ 0 δℓ ν u inequality (34) directly follows, and inequality (33) is now a direct application of the well known formula 2tPp|X| ą tqdt Proof of theorem 3.3. We apply Parseval's formula to derive The second term is easily handled. Remark first that, since s ą 1{2, we have 2s ν`1 ą 2s s`ν`1{2 and we can write ÿ In order to lighten the notations, we will only consider the indexes ℓ ě 2 in the first term. This is of course not problematic, since an identical reasoning allows to bound the two remaining summands by the desired rates of convergence. We hence write the following decomposition ÿ Hence we can decompose further I as Let us first treat the term V I. From Lemmas 6.2 and 6.3 and Cauchy-Schwarz inequality, we derive hich is less than the desired bound for τ large enough. As for term V , we split it in two and write ote ℓ δ "`δ| log δ|˘´2 2ps`νq`1 and write The ε-term is treated similarly by taking ℓ ε "`ε a | log ε|˘´2 2s`2ν`1 and leads to the desired convergence rate. As for the term II, a similar reasonning leads to The term V III is handled as the term V I. Indeed, which is less than the desired rate for τ large enough. Finally, we have Moreover, Lemma 6.1 entails P`B c ℓ˘ď δ κ 2 ρ 2 for all ℓ ě 1. It is hence clear that for κ large enough, the term III is less than the announced rate.
It remains to put together the bounds of the four terms above to get the desired rates of convergence in theorem 3.3.
Proof. The proof is very similar to Lemma 6.3, whence we will just mention the notable changes compared to it. Once more, we shall only treat the case ℓ ě 2. First, we have so that a brief contitionning argument, combined with (30) and Assumption 2.3 entails In order to treat the term Pp|xpK ℓ δ q´11 A ℓ 1 B ℓ δB ℓ f ℓ , ϕ ℓ y| ą tq, we first establish a useful result for the sequel: if g satisfies (20), then }pK ℓ q´1f ℓ } " }T ℓ pf qp1{ 9 gq ℓ } ď }T ℓ pf q} op }p1{ 9 gq ℓ } Furthermore, thanks to Proposition 3.1 we have since f P W s pMq and s ą 1{2. We derive that Let us now bound the term of interest. Once more, we decompose it as r 1`r2 where r 1 and r 2 are defined in (36). We now apply Proposition 16 and (32) and derive The latter is a gaussian random variable with variance δ 2 }pK ℓ q´1f ℓ } 2 À δ 2 ℓ 2ν´1 where we used (40). Turning to the term r 2 , we apply Proposition 16 to derive e now apply (30) and (40) to get }pK ℓ δ q´1f ℓ }} t pK ℓ q´1ϕ ℓ } À ℓ 2ν´1 Hence, et us take a look back to Lemma 6.2. On A ℓ X B ℓ we have prooved that . We deduce À exp`´t δℓ ν´1{2˘1 ttąβ 2 0 δℓ ν´1{2 u`1ttďβ 2 0 δℓ ν´1{2 u The end of the proof is identical to Lemma 6.3.
Proof of Theorem 4.5. The proof is very similar to Theorem 3.3, whence we will just emphasize the notable changes compared to it. First, we apply Parseval's formula to derive The second term is easily handled, since ÿ Thanks to Lemma 6.2 and the definition of S II ℓ,ε , we have on A ℓ X B ℓ . Thus, the Terms I and II can be treated identically to the preceding proof and yield the desired rates of convergence. The terms III and IV are treated exactly as in the preceding proof.

Proof of theorem 4.6
Proof. The lower bound will not decrease for increasing noise levels δ and ε, whence it suffices to provide the case δ " 0 and the case ε " 0 separately. In the sequel, c i will denote a positive constant to be adjusted later and L will play the role of a maximal level. Also, we will note K ν (resp. 9 g ν ) the operator (resp. the function) associated with the Laurent serie p1´zq ν L . The function 9 g´1 {2 will play an essential role in the sequel. Unfortunately it is not square integrable. We thus begin with a preliminary lemma, which states that a minor modification corrects this defect. Lemma 6.5. Let h be the function associated to the Laurent serie Then h is square integrable. Furthermore, for all ν ě 0, for all ℓ ď L, pℓ _ 1q ν´1{2 logpℓ _ 2q À |xK´νh, ϕ ℓ y| À pℓ _ 1q ν´1{2 Proof of Lemma 6.5. h is trivially squared integrable thanks to (26)  This entails the result.
‚ Case δ " 0. For more clarity, we will suppose that ξ is a white noise (the proof readily adapts otherwise). Let hence K 0 " c 1 K ν . Then K 0 P G ν pQq for an appropriate constant c 1 , thanks to Proposition 4.2. Following the arguments of Willer [28], it suffices to find f 0 , f 1 such that i) f 0 , f 1 P W s pMq ii) }f 0´f 1 } 2 Á ε 2s s`ν | log ε|´2 iii) KpP 1 , P 2 q À 1 where P i is the law of y under the hypothesis f i , and K is the Kullback-Leibler divergence.
‚ Case ε " 0. Let L " c 1 δ´1 s`ν . Following the lines of Hoffmann and Reiß [17], we set f 0 " c 2 ϕ 0 , K 0 " c 3 K ν and we only consider couples pK, f q such that Kf " q 0 for a fixed q 0 " K 0 f 0 . It is clear that, for well chosen c 2 and c 3 , we have f 0 P W s pMq and K 0 P G ν pQq. We thus define H the operator associated to the kernel h and introduce K δ " K ν`c4 δH a perturbation of K ν . We shall refer to 9 g δ for the corresponding kernel. Remark that we have f 1´f 0 " c 4 δpK δ q´1Hf 0 " c 4 δpK δ q´1h (41)


extfieldsofstudy:
Computer Science
Mathematics

text:
The straight line complexity of small factorials and primorials

In this paper we determine the straight-line complexity of $n!$ for $n\leq 22$ and give bounds for the complexities up to $n=46$. In the same way we determine the straight-line complexity of the product of the first primes up to $p=31$ and gives bounds for $p\leq 43$. Our results are based on an exhaustive computer search of the short length straight-line programs.

Introduction
In [SS95] Shub and Smale studied the complexity of a number of different algebraic problems in terms of the number of ring operations needed to compute a given ring element by a straight line program. A straight line program for for an integer y can be described as a sequence of tuples x k = (x i • x j ), where i ≤ j < k, x 1 = 1, • can be any of +, −, ×, and the final element x f is equal to y. The smallest integer f such that there exists a straight line program of length f is called the straight line complexity, or cost, of y and is denoted by τ (y) In [BSS89] a general complexity theory for computation over rings was introduced, see also [BCSS98], and the results from [SS95] implies that if the ultimate complexity of n!, here denoted as τ ′ (y), defined as the minimum τ (y) for all y which are integer multiples of n!, is less than (log n) c for some constant c, then P = N P for computation over the complex numbers, and otherwise P = N P . The existence of such a constant c would among other things provide a fast algorithm for factoring integers, see the discussion in [Che04, Che03], and the non-existence would provide strong lower bound for several important problems in complexity theory, see [Bür09] and [Koi04].
The results of [SS95,dMS96,dAM97] provide upper and lower bounds for the straight line complexity of general integers and imply that for most integers τ (n) is not O(p(log log n)) for any polynomial p. In [LSZ82] similar bounds were derived for functions over finite fields. The known bounds for a general integer n are log 2 (log 2 n) + 1 ≤ τ (n) ≤ 2 log 2 n (1) The lower bound is optimal since τ (2 2 k ) = k + 1. The upper bound is achieved by first computing the necessary powers of 2 and then adding them according to the binary expansion of n For specific integers, such as n!, there are few results which strengthen the general bounds. However for n! Cheng derived an improved algorithm, conditional on a conjecture regarding the distribution of smooth integers, and earlier [Str77] a weaker, unconditional, bound was derived by Strassen.
The purpose of this short note is to report the exact values of τ ′ (n!) for small values of n and likewise for τ ′ (p#), where p# is the primorial, which is the product of all primes less than or equal to p. It is easy to see that given a short straight line program for p# we can also find one for n! by using repeated squaring. Our results were obtained by first doing an exhaustive computer search of all straight line program up to a given length, followed by an extended search adapted to finding program for n! and p#.
Most of the material in this note was originally part of a longer paper but while preparing that paper the author found out that the non-computational results were already covered by other recently published papers. That was nearly ten years ago but given the slow progress on problems in this area we hope that these exact results and bounds will help draw attention to the problems and stimulate interest among new researchers. The only additions to the material from the older paper is a recomputation of all data using a newly written program and as a result of this an improvement of some of the lower bounds.
2 Searching for optimal straight line programs  2. Reduce the set of all such programs by removing redundant programs 3. Repeat step 1 and 2 for as large k as possible.
In step 2 we removed all programs in which x k was a number which was also x j for some j < k. For every sequence of numbers (x 1 , . . . , x k ) we kept only one of the programs which computed a permutation of this sequence. We also deleted all program which had x k ≤ 0. The computer resources at hand allowed us to continue this process up to k = 9. In Figure 1 we display some statistics for these programs. We say that an integer y has been reached if there is a straight line program of length at most k which computes y and that y has been covered if y is a divisor of x j for some j ≤ k. We also include the length of the longest interval of the form [1, . . . , x] in which all integers have been reached or covered. After these programs had been found a second search stage was done, with the specific aim. For each n such that τ ′ (n!) ≥ 10 each program of length 9 was extended in a depth first search. During the search every program in which the largest number was not large enough to given a multiple of n! was pruned from the search, as were methods which could not be optimal due to repeated values. Here each program was extended to k = 12, but the number of programs was too large to save so only those which computed a multiple of n! were recorded. Thus all relevant programs of length 12 have been searched, providing us with a lower bound for τ ′ (n!).
The optimal methods are noticeably better than then upper bound for τ (n!) given in 1. The method of Strassen [Str77] gives a bound τ (n!) = O( √ n log 2 n), which seems to deviate more and more from the optimal methods for larger n. The conditional method of Cheng [Che04] has a complexity of the form O(exp(c √ log n log log n)), which certainly seems compatible with the results for small n, but is so sensistive to the value of the constant c that very little can be said based on small values of n.
The same depth first search was done in order to extend the optimal programs of length 12 found for each factorial, as well as few heuristically chosen ones, in an attempt to find short programs for larger factorials as well. In Figure  2 we show the exact values for τ ′ (n!) for n ≤ 22 and for each such n an example of an optimal straight line program. For larger n we display the best method found by our partial search. The final columns states whether the method is optimal or not, and otherwise the lowest possible value.
Finally a similar extension search was performed in order to find optimal methods for the primorials p#. The resulting values and programs are displayed in Figure 3


text:
Sparse Penalty in Deep Belief Networks: Using the Mixed Norm Constraint

Deep Belief Networks (DBN) have been successfully applied on popular machine learning tasks. Specifically, when applied on hand-written digit recognition, DBNs have achieved approximate accuracy rates of 98.8%. In an effort to optimize the data representation achieved by the DBN and maximize their descriptive power, recent advances have focused on inducing sparse constraints at each layer of the DBN. In this paper we present a theoretical approach for sparse constraints in the DBN using the mixed norm for both non-overlapping and overlapping groups. We explore how these constraints affect the classification accuracy for digit recognition in three different datasets (MNIST, USPS, RIMES) and provide initial estimations of their usefulness by altering different parameters such as the group size and overlap percentage.

Introduction
Restricted Boltzmann Machines (RBMs) are Energy Based Models (EBMs) that have been extensively used for a diverse set of machine learning applications mainly due to their generative and unsupervised learning framework. These applications range from image scene recognition and generation [? ], video-sequence recognition [? ] and dimensionality reduction [? ].
An equally important aspect of RBMs is that they serve as the building blocks of DBNs [? ]. Their use as such has been favored in the machine learning community due to the conditional independence between the hidden units in the RBM that allows for the efficient and computationally tractable implementation of deep architectures.
In recent years, sparsity has become an important requirement in both shallow [add cite] and deep architectures. Although primarily used in statistics for optimization tasks in order to overcome the curse of dimensionality in various applications, it also serves as a way to emulate biologically plausible models of the human visual cortex, where it has been shown that sparsity is an integral process in the hierarchical processing of visual information [? ? ? ].
Moreover, an added benefit of using sparse constraints in the form of mixed norm regularizers in deep architectures is that they can alleviate their restrictive nature by allowing implicit interactions between the hidden units in the RBMs. Mixed norm regularizers such as l 1,2 have been extensively used in statistics and machine learning [? ]. In this paper we provide initial results when inducing sparse constraints by using a mixed norm regularizer on the activation probabilities of the RBMs. The mixed norm is applied on both non-overlapping and overlapping groups. We also show that this regularizer can be used to train DBNs, and offer results for the task of digit recognition using several datasets.

Restricted Boltzmann Machines
An RBM is a type of two layer neural network comprised of a visible layer that represents the observed data x and a hidden layer that represents the hidden variables h. The addition of these hidden units allows the model an increased capacity in expressing the underlying distribution of the observed data.
RBMs are energy based models and as such they define a probability distribution through an energy function as seen in Eq. 1 Where Z, provided in Eq. 2, is called the partition function and is a normalizing factor ensuring that Eq. 1 is a probability.
In the case of an RBM the energy function E(x, h) is defined in Eq. 3: b are the visible unit biases and a are the hidden unit biases.
In the common case where we are using stochastic binary units for both visible and hidden units, then the conditional probabilities of activation are obtained by: where σ is the sigmoid function and Since an RBM does not allow for connections amongst hidden units or amongst visible units we can easily obtain Eq. 6.
Intuitively, the observed data, x will be modeled by those hidden units, h that are expressed with a high conditional probability p(h j |x). The goal of adding sparse constraints to the network is to allow for the salient activation of the hidden units based on the differences of the observed data. As a result, we can achieve an initial clustering of the observed data that will increase the discriminative power of the model.

Training an RBM
RBMs are energy based, generative models that are trained to model the marginal probability p(x) of the observed data where: In general, energy based models can be learnt by performing gradient descent on the negative loglikelihood of the observed data. Specifically, to learn the parameters of the network we need to compute the gradient provided in Eq. 8 given the observed (training) data x l .
where · n denotes the expectation with respect to n. As evident in Eq.8, the gradient has two phases. The positive phase which tries to lower the energy of the training data x l and the negative phase which tries to increase the energy of all x in the model. Assessing the energy on all the data can be an intractable task given the size of the network and the number of possible configurations. In order to obtain an approximation Hinton (2006) successfully proposed the use of Contrastive Divergence (CD). This allows us to sample an approximation of the expectation over (x, h) using Gibbs sampling at only k steps. Empirically, it has been shown that setting k = 1 will provide an adequate approximation although it will not follow the theoretical gradient [? ].
Applying CD on Eq. 8 we can obtain the following update equations for the parameters of the network.
where the ( ·) defines the generated distributions obtained by the CD.
In the next section, we introduce a general version of sparse constraints in the learning phase of the RBM through the use of the mixed norm in an effort to control the activation probabilities of the hidden units.

Mixed Norm RBMs
Several attempts in inducing sparse constraints in the RBM by [? ? ] have been successful in increasing the discriminative power of the models. Examples of these sparse constraints range from weight decay [? ] to modified norm penalties [? ]. In this paper we focus on the generalized penalty of the mixed norm (l 1,2 ), but will also provide a theoretical and practical implementation for the use of overlapping groups. We will refer to this generalized penalty applied to the expectations of the activation probabilities as the Mixed Norm RBM (MNRBM).
As mentioned before, learning an RBM consists of performing gradient descent on the negative log-likelihood. We can thus define the cost function L to be minimized as L = −logp(x). When applying the mixed norm regularizer the cost function takes the general form of Eq. 12.
Where λ is a regularizer constant. The second term of Eq. 12 defines the mixed norm penalty on the expectations of the hidden unit activation probabilities. In order to apply the mixed norm we assume that the hidden units are divided into groups. These groups can be non-overlapping or overlapping. As a result, we are able to penalize a whole group and not just individual hidden units.
MNRBM with non-overlapping groups: Given an RBM with J hidden units we define a partition of the hidden units into groups P m where m = 1, 2, ...M . The groups are non-overlapping and of equal size to alleviate computational issues. The mixed norm penalty for a data sample x l is defined in Eq. 14.
In practice, the desire behind the application of the mixed norm penalty is to set groups of the hidden units to zero when representing the observed data by forcing their activation probabilities to zero. As a result, given an observed data sample only a small number of groups of hidden units will be activated, leading to its sparse representation.
MNRBM with overlapping groups: Given an RBM with J hidden units we define a partition of the hidden units into groups P m where m = 1, 2, ...M . The groups are overlapping and of equal size. Depending on the percentage of overlap, a we will obtain a new set of groups P k where k = 1, 2, ...K.
We can then define a set of augmented hidden units J = {h ∈ J : ∀P k , P k P m = P m , J ⊃ J}. Subsequently, given an RBM with J hidden units, we can then consider that the set P k , defines non-overlapping, equally sized groups [? ]. The mixed norm penalty for a data sample x l is defined in a similar way as in Eq. 14.

Training the Mixed Norm RBM
In order to train the MNRBM with non-overlapping groups and obtain the model parameters θ we need to minimize the cost function presented in Eq. 12. This can be achieved by performing a coordinate descent once we have obtained the gradients of the regularizers.
The gradient of the mixed norm penalty for the weights, W is as follows: When applied on the expectations of the activation probabilities the mixed norm penalty will follow their trend while forcing the groups that include members with low activation probabilities towards zero. The l 2 norm in the denominator ensures that the groups with low activations will be pushed further closer to zero.
Given the gradients of the penalties the update equations for the MNRBM are presented bellow: The detailed steps for training the MNRBM are depicted in Algorithm 1.
The general penalty of Eq. 12 allows us through the manipulation of the constant regularizer, λ, the group size and percentage of overlap to obtain different types of architectures. In this case, the sparsity is induced at the group level of the hidden units whereby the observed data is represented by a small number of groups of hidden units. The λ constant is empirically determined [? ] based on the task at hand. However, there appears to be a bimodality whereby a large proportion of the activation probabilities is set to a high value for the non-overlapping groups MNRBM. This may be attributed to the choice and size of groups when applying the mixed norm penalty. Given that the activation probabilities are pushed towards high values one can expect that such a process may have an adverse result for classification tasks since the hidden units will over-represent the observed data. However, in the case of the overlapping groups most of the activations are pushed towards zero. Although, that is the goal of adding sparse constraints is to force the activations to zero, in this case we may be dealing with a biased system that actually under represents the data distribution.

Data
We have used three different data sets in order to train and test the network.
• MNIST is a popular data set in the community for hand-written digit recognition and is comprised of 70000, 28 × 28 images (60000 train -10000 test). It is publicly available at yann.lecun.com/exdb/mnist.
• The RIMES data set which was created by asking volunteers to write hand written letters for different scenarios. In this paper we used the digit set of the data base. In total the set we used was comprised of 37200 images of different sizes (29800 train -7400 test). Further information can be obtained at www.rimes-database.fr.
• The USPS digit data set that we used is comprised of 9280 (7280 train -2000 test), 16 × 16 images. The extracted images were scanned from mail in working U.S. Post Offices [? ].
In order to achieve the cross-training and testing all images were resized to have the same size as the MNIST dataset (28 × 28) given its extensive use in this task. All images were also checked to ensure that orientations/translations were uniform across the data sets. No other pre-processing was employed. Example images from the three datasets can be seen in figure 3 In order to offer a comparative view between the different architectures we used Hinton's model for digit recognition, but we substituted the vanilla RBM with the proposed MNRBM. We pre-trained a 500 − 500 − 2000 DBN and tested it on three different data sets, MNIST, RIMES and USPS.
Continuing, to obtain classification error rates we added 10 softmax layers to get the posterior probabilities for the different classes. The network was fine-tuned using conjugate gradient as described in [? ]. The constant regularizer was empirically set to λ = 0.1 for all different models [? ].
Continuing, for the mixed norm architecture with non-overlapping groups we used different group sizes for the hidden units, 5, 10, 20 and 100 respectively. In the case of overlapping groups we used group sizes of 20, and 50 with a = {20%, 50%}. Results on the classification accuracy and the computational cost of the models can be seen in Table 1 and Table 2 respectively. All experiments were performed on a 24 core server (AMD Opteron processor 8435) with a core CPU of 2.6GHz and a cache of 512KB.
(a) Average Probability density function for activation probabilities, p(h = 1|x) (b) Classification Accuracy for the different penalties using the USPS data set Figure 4: Average log pdf for activation probabilities for the vanilla RBM and mixed norm RBM using a batch of the USPS data set and classification accuracy for the USPS data set using the different architectures From Table 1 we can infer that our proposed mixed norm penalty can offer the flexibility of creating architectures that will be able to matchthe classification accuracy of the models depending on the underlying distributions. It appears that for the task of hand-written digit recognition the distribution of the observed data favors the use of larger non-overlapping group sizes for the mixed norm architectures.
In order to get a better understanding of the impact of the different sparse constraints and architectures, Figure 4 depicts the average probability density functions of the expectations of the activation probabilities for the MNIST training data.
It is interesting to note that the proposed architectures that utilize the mixed norm penalty (MNDBN)with the overlapping groups tend to aggressively push their activation probabilities to zero. However, these architectures also tend to offer lower accuracy rates which can be attributed to an inability of the models to concisely capture the underlying data. A possible way for exploring this phenomenon further may be to constrain the penalty of the expectations as seen in [? ].

Conclusions
In this work we provided some first insights for the use of the mixed norm sparse constraint in DBNs. We performed experiments using three different data sets for the task of hand written digit recognition and offered a practical approach for the use of overlapping groups with the mixed norm constraint. Although, our initial experiments were limited in the use of equal size overlapping groups, one could easily extended to non-symmetric overlapping groups, using a similar methodology. Inducing sparse constraints based on specific geometries may also provide better results in the case of digit recognition and offer more interesting results for tasks such as scene categorization.


extfieldsofstudy: Psychology Medicine

text:
An improved algorithm for model-based analysis of evoked skin conductance responses☆

Highlights • We improve predictive validity of a general linear convolution method to analyse evoked SCR.• A constrained individual response function provides highest predictive validity.• This IRF is realised by a canonical SCRF together with its time derivative.• A high pass filter of 0.05 Hz cut-off frequency is optimal for analysis.• Non-linear models better reconstruct the observed time-series but have lower predictive validity.

Introduction
Recent interest in model-based analysis of skin conductance responses (SCR) (Bach & Friston, 2013) is -in part -motivated by the need to improve the temporal resolution of inference in rapid event-related paradigms (Barry, Feldmann, Gordon, Cocker, & Rennie, 1993). In model-based approaches, generative (forward) models specify how underlying physiological or psychological states generate observed data. Model inversion refers to estimating these (hidden) states from data. It turns out that inversion of probabilistic forward models has fundamental advantages, one of them being a propensity to suppress the effect of measurement noise (Bach, Flandin, Friston, & Dolan, 2009;Bach, Daunizeau, Friston, & Dolan, 2010;Bach, Daunizeau, Kuelzow, Friston, & Dolan, 2011). Statistical inference on the hidden states is generally more powerful than statistical comparisons of observed data because the models are more informed or constrained, leaving greater degrees of freedom in the data for efficient inference. Furthermore, the parameters of generative models provide a quantitative and explicit description of assumptions implicit in operational approaches (Bach & Friston, 2013), thus allowing for rigorous testing of those assumptions. Finally, model-based approaches afford quantitative rather than qualitative description of hidden, psychological processes.
Evoked skin conductance responses (eSCR) that follow a short (less than second) stimulus can be analysed with general linear convolution models -similar to the convolution models widely used in the analysis of functional magnetic resonance images (Friston, Jezzard, & Turner, 1994). In order to estimate the amplitude of evoked sympathetic nerve activity (SNA) from eSCR, we proposed such a convolution model (Bach et al., 2009). This model comprises two parts: a peripheral model incorporating a (standard linear time invariant) convolution operator, thought to be implemented by the sudomotor nerve terminal and sweat gland; and a linear neural model assuming infinitely short neural bursts immediately after each stimulus. We have shown that time invariance assumptions for the peripheral system are largely met (Bach, Flandin, Friston, & Dolan, 2010), while non-linearities in the peripheral system may occur but can easily be modelled within this framework -see (Bach & Friston, 2013) for a discussion. The model is highly regularised by placing informative constraints on the form or shape of the convolution kernel which models the peripheral response function (RF). This enables one to estimate the mean evoked SNA amplitude for each condition of an experimental design -even when observed eSCR overlap in time. This model was designed to optimise the predictive power of the estimates, rather than to precisely reconstruct the observed time series. Indeed, when subjects observe negative-arousing or neutral pictures, picture category can be better predicted from SNA estimates than from SCR peaks, an observation that speaks to its predictive validity (Bach et al., 2009).
Model-based eSCR analysis, based on probabilistic inversion of a general linear convolution models, is thus a potentially powerful method. As with any method, however, the practical implementation makes certain technical assumptions that go beyond the known biophysical properties of the system. Three points deserve particular attention: Firstly, the peripheral response model uses a canonical skin conductance response function (SCRF) for all experiments and individuals. Such a stereotypical response function is a strong biophysical assumption and unsupported by observation. Indeed in our own validation experiments, we observed large interindividual variability, accounting for up to 20% of overall response variability (Bach, Flandin, et al., 2010). Therefore, we added orthogonalised Taylor expansions to the SCRF to account for differences between individuals and conditions, thus improving model fit (Bach et al., 2009). Effectively, this enables the model to fit a subject specific RF in terms of a linear mixture of basis functions of peristimulus time, where the basis set is generated by the Taylor expansion. However, because the additional basis functions are orthogonalised to the SCRF, they do not affect the estimation of the parameter controlling the amplitude of the SCRF (Calhoun, Stevens, Pearlson, & Kiehl, 2004;Hopfinger, Buchel, Holmes, & Friston, 2000). Yet, this parameter is taken to estimate the SNA. This means that additional basis functions improve data fit at the within subject level but not comparisons of SNA at the between subject level. Hence, one might ask whether modelling an individual response function (IRF) -rather than a canonical skin conductance response function (SCRF) -improves predictive validity.
There are several ways to model subject specific IRF. First, the SCRF together with the remaining basis functions can be used to estimate a subject and condition specific IRF. That is, we can reconstruct the estimated eSCR, measure the peak amplitude (over peristimulus time) and use this as an estimate of SNA amplitude, instead of the canonical parameter estimate. Other regularised basis sets also provide models of IRFs. An uninformed finite impulse response (FIR) model was proposed in Bach et al. (2009) due to its popularity in fMRI research. A cosine set also used in fMRI analysis serves the same purpose. These basis sets typically have a larger number of basis functions than basis sets built upon truncated Taylor expansions. This means that although they are more flexible, they require greater numbers of parameters to be estimated. In all these approaches, separate IRFs are estimated for each condition within one individual. A more informed approach is to assume the form of the subject specific response function is the same for all experimental conditions. We have implemented this constraint by extracting data from all conditions and fitting a response function to the first principal component of the data. We will refer to this response function as the subject-specific response function (SRF).
A second issue is that skin conductance time series comprise both phasic responses and a slowly drifting tonic component. This drift is why many analysis schemes high pass filter the signal (Boucsein, 2012), including contemporary model-based approaches (Benedek & Kaernbach, 2010a, 2010b. This renders the phasic responses finite and removes slow signal drifts which are difficult to model. In our implementation, we used a bidirectional firstorder Butterworth filter with time constant of 10 s (corresponding to a cut off frequency of 0.0159 Hz) (Boucsein, 2012). A bidirectional filter was chosen as it retains peak latencies. Because this filter can slightly distort the shape of the signal, the regressors of the general linear convolution model are subjected to the same filter. The choice of the filter frequency is based on prior experience but not on theoretical considerations. Therefore, one may ask whether there is an optimal filter that provides the best data features for modelling. Generally speaking -when modelling biological time series -data features that cannot be produced by a plausible forward model are probably measurement noise or the product of hidden processes not included explicitly in the model. This usually means they can be discarded with impunity, thereby increasing the signal-to-noise ratio (SNR) of the pre-processed data. Data conditioning is then, effectively, a part of model inversion. The question here is whether there is an optimal high pass filtering of skin conductance timeseries that increases signal-to-noise. In case of a signal with precisely known RF, the matched filter theorem provides a way of theoretically deriving a filter that maximises the SNR. In our case, the true RF is not precisely known, and also varies between individuals, such that we sought to empirically determine the filter characteristics that maximise predictive validity of SN estimates.
Finally, a linear neuronal model makes the strong assumption that SNA evoked by a short stimulus occurs with constant latency. We have previously shown that under this assumption there is no evidence for time-varying responses in the peripheral system (Bach, Flandin, et al., 2010). Here, we revisit this assumption and investigate whether modelling variations in neuronal latency improves predictive validity, under the assumption of an invariant peripheral response. Hence, we compare linear and non-linear models. Two particular non-linear models are considered. First, we used our previous approach that uses Dynamic Causal Modelling (DCM) (Bach, Daunizeau, et al., 2010) to obtain estimates of SNA amplitude per trial by letting response amplitude and onset vary on a trial-by-trial basis. Note that the neural model here is still informed insofar as it specifies a certain response window. Some authors propose uninformed neural models; in other words, they assume that SNA can occur any time, but to only use SNA during post-stimulus time windows for analysis (Benedek & Kaernbach, 2010a, 2010b. We sought to emulate this approach using DCM for spontaneous fluctuations (Bach et al., 2011). Both approaches yield a trial-by-trial estimate of SNA amplitude, which was averaged across experimental conditions for comparison with other approaches.
In some circumstances -e.g. to use neural response estimates as explanatory variables for analysis of independent experimental data, such as fMRI, trial-by-trial estimates of SN amplitudes may be required. Here, we sought to establish whether linear models are sufficient for this purpose or whether the iterative procedures required for non-linear models inherent in DCM are justified.
We have previously discussed how to benchmark methods that estimate hidden variables from observed data (Bach & Friston, 2013). One way is by making certain assumptions about what causes the hidden variable to change. A consensus assumption in the psychophysiology literature is that emotionally arousing events increase sympathetic arousal, as engendered by negative and positive arousing images. This has been demonstrated using operational approaches (Amrhein, Muhlberger, Pauli, & Wiedemann, 2004;Greenwald & Lang, 1989;Johnsen, Thayer, & Hugdahl, 1995;Winton, Putnam, & Krauss, 1984). Here, we assume that negative and positive arousing images would elicit greater sympathetic arousal than neutral non-arousing images, and evaluated different models in terms of their ability to distinguish between image types, using just the observed SCR.
In summary, we evaluated our method empirically, by comparing the predictive validity of different generative models. In a first step, we compared a canonical response function against various forms of an individualised response functions (IRF). Taking the best model from this step, we then compared various filter settings, non-linear methods, and the efficiency of trial-by-trial estimates. These comparisons used two independent, previously unpublished data sets.

Datasets
We analysed two datasets, both acquired during a study that examined the influence of distracting background noise on the perception of emotionally arousing images. There was no interaction of acoustic distractors with differential SCR in either experiment. The main effects of auditory stimulation will be reported elsewhere. For the present analysis, all data were collapsed across the auditory stimulation factor. Participants for both experiments were recruited from UCL students and in the general population via an online recruitment system. Both experiments were approved by the Joint UCL/UCLH Committees on the Ethics of Human Research.
The first experiment investigated SCR in response to negative-arousing and neutral pictures. 60 healthy individuals (30 male, 30 female; age: M = 23.7; SD = 4.7 years) participated in the first experiment. Participants watched, in randomised order, the 45 least arousing neutral (valence within 1 standard deviation around the mean) and 45 most arousing aversive pictures (valence lower than 1 standard deviation below the mean) from the International Affective Picture Set [IAPS] (Lang, Bradley, & Cuthbert, 2005) for 1 s each, with an inter stimulus (ISI) interval randomly determined as 7.65 s, 9 s, or 10.35 s. Participants were instructed to press the cursor up or down key on a computer keyboard to indicate whether they liked the picture or not. This response served to increase SCR. The experiment was divided into three blocks with 45 s breaks in between.
The second experiment was similar to the first, with an additional third condition of positive arousing images. 40 healthy individuals (20 male, 20 female; age: M = 21.9; SD = 3.8 years) watched the 16 least arousing neutral, most arousing aversive and most arousing positive images (defined analogous to experiment 1, and excluding explicit nude images) from the IAPS, for 1 s each, in randomised order, and in one single block. ISI was 4.4 s. Responses were the same as in experiment 1.
Lists of images, used in both experiments, are available from the authors.

SCR recordings and preprocessing
We recorded skin conductance on thenar/hypothenar surface of the nondominant hand using 8 mm Ag/AgCl cup electrodes (EL258, Biopac Systems, Goleta CA, USA) and 0.5%-NaCl electrode paste (GEL101; Biopac Systems), using a custombuild constant voltage coupler (2.5 V). The output of the coupler was converted into an optical pulse frequency with an offset (i.e. minimum sampling rate) of 100 Hz and a factor of 31 Hz/S; the converter was designed for a maximum conductance of 50 S, corresponding to 1650 Hz pulse frequency. The optical signal was converted to voltage pulses and recorded (Micro1401/Spike 2, Cambridge Electronic Design, Cambridge, UK). SCR data were filtered with a 1st order Butterworth filter. A bidirectional filter with low pass cut off frequency of 5 Hz was used, and a bidirectional high pass filter with cut off frequency of 0.0159 Hz for step 1. High pass filter settings were variable for steps 2-4. Data were then down sampled to 10 Hz.

Inversion schemes and models
In steps 1 and 2, the data were analysed with a general linear convolution model as proposed in (Bach et al., 2009). Additional models were compared in step 3-4, see below. For each experimental condition, event onsets were modelled as Dirac delta functions and convolved with a single RF or basis set as described below. The resulting design matrix of the general linear model contained one regressor for each condition and function in the basis set. For basis sets containing more than one basis function, the amplitude of the underlying SNA was estimated as the peak of largest modulus, of the reconstructed SCR from all basis functions in the basis set. To ensure that these findings could be generalised from the detection of categorical conditionspecific differences to continuous variables, we additionally analysed a model for experiment 1 in which picture onset was modelled by a single stimulus function for both conditions, and picture type and linear habituation effects were modelled as parametric modulators. This models the linear component of the (possibly nonlinear) habituation of SNA amplitude over the course of an experiment.

Step 1: response functions
The basis sets we assessed are listed below. The number of basis functions (bf) and of free parameters (k) for a model containing j conditions and N trials are noted in brackets for each basis set. For basis sets containing more than one basis function, the ensuing regressors were orthonogalised using a serial Gram-Schmidt procedure: (0) Peak scoring: several peak scoring methods were computed to allow comparison with the literature. We followed the recommendations of the Society for Psychophysiological Research (SPR)  to identify an SCR according to the onset latency, identified by the point of maximum deflection, together with the rise-time (i.e. onset-peak latency) of the subsequent peak. Because there is no community consensus on the optimal duration of the onset latency window, we used a window of 1-4 s (Boucsein, 2012;Dawson & Filion, 2007;Edelberg, 1972) and a window of 1-3 s (Barry, 1987;Boucsein, 2012;Dawson & Filion, 2007), and a post-onset peak window of 0.5-5 s (Boucsein, 2012). Onset SCR values were subtracted from peak SCR values. Non-responses were scored zero. For both windows, we calculated SCR amplitude by omitting responses below 0.01 S before averaging, and SCR magnitude by averaging all responses including zero responses. Further, we used a simpler algorithm and subtracted the mean of a 1 s pre-stimulus from the highest value in a 1-4 s post-stimulus response window.
(1) SCRF: the standard SCRF (bf = 1, k = j) -this is the reference model for steps 1 and 3. (2) SCRF/time derivative: the standard SCRF with time derivative (bf = 2, k = 2j)this is the reference model for step 2. to the first principal component with a 2nd order ordinary differential equation, using a variational Bayes inversion scheme as described in (Bach, Daunizeau, et al., 2010). This form of model was considered for experiment 1 -in experiment 2, the minimum SOA was too short to estimate a robust SRF (bf = 1, k = 3 + j).

Step 2: filter settings
In the second step, we compared different filter settings for the best method from step 1 (SCRF with time derivatives) which we define as reference method for step 2. Uni-and bidirectional 1st order Butterworth high pass filters were compared with cut-off frequencies of 0.005 Hz, 0.1 Hz, 0.0159 Hz (the current default), and from 0.02 Hz to 0.10 Hz in 0.005 Hz steps. Bidirectional filtering, as implemented in the Matlab function filtfilt.m, filters the time series twice, in the forward and backward direction. To keep filter order precisely the same, we filtered twice for unidirectional filtering, but both times in forward direction. To ensure the robustness of these findings, the evaluation of different filtering was also performed using a simple SCRF without derivatives.

2.6.
Step 3: non-linear models (1) Informed DCM: In the first non-linear model, an SN burst was assumed to occur within 2000 ms after stimulus onset -this differs from the linear model where it was assumed to occur immediately after stimulus onset. The onset, duration and amplitude of each SN burst was estimated using DCM as described in (Bach, Daunizeau, et al., 2010) (k = 3N). (2) Uninformed DCM: In a second non-linear model, SN bursts were assumed to occur anywhere within the experiment (with fixed dispersion to constrain the number of free parameters) -analogous to our DCM for spontaneous fluctuations (Bach et al., 2011). Due to the large numbers of parameters, the model could not be inverted, so a 60 s epoch was inverted every 30 s. From the overlapping epochs, the middle 30 s were analysed. Each estimated SN response that caused an SCR that fell into a 1-4 s post stimulus time window was extracted, and the largest response in each window was retained. This model is entirely uninformed about the experiment, in line with other contemporary modelbased approaches (Benedek & Kaernbach, 2010a, 2010b) (k = 60e, where e is the number of inverted data epochs).
Both non-linear models were compared to the reference model from step 1.

2.7.
Step 4: trial-by-trial estimates All analyses described thus far estimated SNA amplitude under the assumption that the responses were the same for each trial in a condition. In some circumstances one might want to estimate trial by trial SNA. To evaluate the validity of trial by trial estimates, we assumed that arousal ratings from the validation sample of the IAPS (Lang et al., 2005) and habituation were the principal causes of trial by trial variations. We therefore quantified the proportion of variance in trial-by-trial estimates explained by these variables. We analysed single trials from experiment 1 and computed a GLM as explained above, with one regressor per trial, using either the SCRF or SCRF with time derivative, and a high-pass filter with cut off frequency of 0.05 Hz. For comparison, a DCM was inverted where SCR were modelled as evoked responses with constant latency. This is a linear neural model but using a non-linear, iterative inversion scheme. Finally, we used trial-by-trial SNA from the informed DCM from step 3. These four models were compared against each other.

Model comparison
The different models and data features were compared in terms of their predictive validity; i.e., their ability to predict a stimulus class from estimated SNA, for a particular experimental contrast. The contrast of interest for experiment 1 was the difference between neutral and negative-arousing pictures. For experiment 2, we used the contrasts neutral vs. negative-arousing and neutral vs. positive-arousing. To assess predictive validity, we used a general linear model with the contrast of interest as the response variable, and the estimated SN activity as predictor. The design matrix included subject effects. The residual sum of squares RSS was converted to a negative log likelihood value LL, such that smaller LL values indicate a higher predictive validity using the following relation where n is the number of observations. This disregards model complexity, which was the same for all analyses. We report log evidence differences or log Bayes Factors (LBF) -the difference in log likelihood between each model and a reference model, for each step in model comparison. Negative LBF values indicate a model fit that is better than the reference model. An LBF difference larger than 3 is often considered decisive as it corresponds to a p-value of 0.05. We further compared how well the models fit the observed SCR time series given the optimised parameters. This was done for the sake of completeness, although this measure of model performance should be interpreted with caution because (within subject) measures of accuracy, such as this, do not reflect predictive validity at between subject level. We converted the residual sum of squares RSS from the fitted time series into a negative log likelihood value LL using Eq. (1) which was then corrected for complexity to give the Akaike information criterion (AIC) by the relation where k is the number of free parameters in the model. AIC is an approximation to Bayesian model evidence (Penny et al., 2010). For step 4, we converted the relative residual variance into negative log likelihood by Eq. (1), and summed up the likelihood terms across participants. For model fit, we report Log Bayes Factors (LBF)the difference in AIC between each model and a reference model.

Results
The stimuli had a clear effect on estimated SNA as demonstrated by all analyses. For the reference (standard SCRF) method, negative-arousing images elicited higher SNA than neutral images in experiment 1 (t(59) = 4.70; p < .00001), and in experiment 2 (t(38) = 4.68; p < .00001). The difference in activity between positive and neutral images (t(38) = 1.44; p = 0.16) failed to reach significance in the reference method but was significant in other model-based methods.

Step 1: optimising the response function (RF)
We compared all models against our reference model with SCRF alone, and added a peak scoring approach for comparison with the literature. Fig. 1a shows that, for all contrasts and among linear models, SNA estimates from a model including the SCRF and its time derivative best predicted the experimental manipulation; i.e., they had the lowest Log Bayes Factor (LBF). The LBF difference between this and the reference method was decisive for all three contrasts. Estimates from linear models with more degrees of freedom had equal or lower predictive validity (i.e. higher LBF) than the reference. At the same time, the best linear model was better than any of the peak scoring methods.
In order to assess whether this result was due to modelling of dichotomous contrasts between experimental conditions or whether it would generalise to continuous predictors, we modelled condition and habituation in experiment 1 as parametric modulators. Again, an SCRF with time derivative model had best predictive validity, both for the effect condition and the continuous predictor, the linear habituation term (data not shown).
In terms of within-subject model fit, we found a different picture (Fig. 1). More complex models had a massively improved model fit, despite penalising model complexity, while the predictive validity of these methods tended to be worse. This might indicate that the more complex models over-fit noise variance not induced by the experimental manipulation.
In summary, we conclude that a linear model using SCRF and time derivative provides the approach with the greatest predictive validity (among the ones considered here). This appears to be the case over the two experiments with different experimental contrasts. Hence, we used this model for the next evaluation step.

3.2.
Step 2: optimising the filter All filter settings were referenced against the winning method from step 1, which estimates SN amplitude from a reconstructed response function, using a general linear convolution model with a canonical SCRF and time derivative, at standard high pass filter settings -a bidirectional first order Butterworth filter with cutoff frequency of 0.0159 Hz. In Fig. 2, this model hence has a LBF of zero. We varied filter cut-off frequency from 0.005 to 0.1 Hz, and employed a unidirectional and bidirectional filter. Fig. 2 shows that increasingly higher filter frequencies produced SN response estimates with better predictive validity, up to about 0.05-0.06 Hz for the contrast negative > neutral in both experiments, and up to about 0.035 Hz for the contrast positive > neutral from experiment 2. Beyond this frequency, predictive validity did not change much, but tended to get slightly smaller towards higher frequency cut-offs. There was no clear advantage for either uni-or bidirectional filtering. Unidirectional filters tended to be slightly better for higher frequency and vice versa for lower frequency. Note that both uni-and bidirectional filtering applied the filter twice, but in different directions. The same pattern of results was found when the SCRF was used without derivatives (data not shown), indicating that the influence of filter is not due to peculiarities of the response function.
As in step 2, the effect of the filter upon predictive validity was not precisely reflected in first-level model fit. Model fit asymptotically increased with higher filter cut-off frequencies -as one might expect when degrees of freedom are removed from the data (making it easier to fit).
To summarise, higher filter cut-off frequencies produced better predictive validity. However, across the different contrasts, the precise choice of filter frequency had no decisive and unambiguous effect on predictive validity between a frequency of 0.035 and 0.1 Hz.

Step 3
In the third step, we compared linear and non-linear models which model more degrees of freedom in underlying SNA latency. The DCM for event-related responses allows variability in the neural response timing on a trial-by-trial basis and might thus be able to estimate SN amplitude even if neural latencies fluctuate. DCM for spontaneous fluctuations makes even less assumptions about the neural response and can capture different shapes of neural responses, e.g. mono-or biphasic. Fig. 3 shows that the DCM for event-related responses -labelled informed because it is informed about the experimental timinghad decisively lower predictive validity than any linear model for all three contrasts. The uninformed DCM had lower predictive validity than linear models for the contrast neutral > negative, across both experiments. For the contrast positive > neutral, it had approximately the same predictive validity as the best linear model. As in the previous steps, a better model fit under more complex models with more degrees of freedom was not reflected in predictive validity. Step 1: comparison of a linear model for evoked SCR, using different response functions as explained in Section 2. Lower Log Bayes Factors (LBF) indicate higher model evidence for the target model. Upper panel: predictive validity; i.e., ability of estimated SN amplitudes to predict a known sympathetic state, for three contrasts from two experiments, expressed in LBF as negative log likelihood difference between the model in question and a reference model. Several peak scoring methods (PS) are added for illustrative purposes as null models (left of the dashed line). Lower panel: Model Log evidence of the within-subject model, expressed as difference in AIC between the target model and our benchmark model, summed over participants. Abbreviations: PS: SPR (1-4 s) amp -peak scoring amplitude according to the SPR recommendations, using a 1-4 s post-stimulus onset window; PS: SPR (1-4 s) mag -peak scoring magnitude according to the SPR recommendations, using a 1-4 s post-stimulus onset window; PS: SPR (1-3 s) amp -peak scoring amplitude according to the SPR recommendations, using a 1-3 s post-stimulus onset window; PS: SPR (1-4 s) mag -peak scoring magnitude according to the SPR recommendations, using a 1-4 s post-stimulus onset window; PS: peak/baseline -peak scoring magnitude, substracting a 1 s pre-stimulus baseline from the maximum value within a 1-4 s post-stimulus window; SCRF -skin conductance response function (benchmark method); SCRF/time deriv. -skin conductance response function with time derivative; SCRF/time and disp deriv. -skin conductance response function with time and dispersion derivative; FIR 15 s -uninformed finite impulse response function with 15 timebins of 1 s duration; FIR 30 s -uninformed finite impulse response function with 30 timebins of 1 s duration; Cosine nth order -cosine basis set of nth order; SRF -subject-specific response function.

Step 4
The single-trial GLM (with two different RFs) and a DCM for evoked responses performed almost equally well: the mean explained variance proportion across participants was (mean ± standard deviation) 0.15 ± 0.13, 0.16 ± 0.14, 0.17 ± 0.13, respectively, while the non-linear DCM for event-related responses performed slightly worse, with an explained variance proportion of 0.12 ± 0.09. Model evidence was decisively higher for the iterative DCM for evoked responses than for all other approaches. Step 2: comparison of different unidirectional and bidirectional high pass filters, applied to the data before model inversion. Lower Log Bayes Factors indicate higher target model evidence. All Log Bayes Factors are with respect to our current standard filter, a bidirectional first order Butterworth filter with cut off frequency of 0.0159 Hz.

Discussion
Model-based methods for analysis of SCR responses estimate underlying SNA and are potentially important for non-invasive measurements of autonomic responses. Here, we sought to validate and improve a method based upon linear and nonlinear convolution models -with three results.
First, as general validation of linear models, the best linear method had as good, or better, predictive validity than various peak-scoring approaches based directly on observed data, or than all non-linear methods. At the same time, first-level model fit was unambiguously higher for models with higher complexity, in particular for non-linear models. This demonstrates that the tight constraints imposed by linear models preclude overfitting, and to help recover the underlying hidden cause more efficiently.
Second, the original linear method was improved by allowing for between-subject variability in the form of peripheral responses. The best predictive validity was achieved with a model that used a standard canonical SCRF across subjects, together with its temporal derivative, to estimate an individual and condition specific RF. It is noteworthy that the approach was also best when the binary variable (experimental condition) was encoded as parametric modulator, or when predicting a continuous variable (linear response habituation) encoded as parametric modulator. This approach of complementing the SCRF with Taylor expansions was already implicit in the original model proposed in (Bach et al., 2009;Bach, Flandin, et al., 2010); but in the original method it was only used to improve within-subject model fit, not to estimate SNA amplitude. Here, we extend this approach and propose this as a future standard for analysis of eSCR. Note that this remains a constrained or informed way of estimating IRFs. Other approaches with even slightly higher complexity did not improve predictive validity. In a nutshell, the winning method uses the same approach that has been proven most sensitive in fMRI analysis (Hopfinger et al., 2000). It is interesting to note that a more constrained approach, namely the estimation of a subject-specific response function (SRF) from the data without allowing for condition differences, proved less powerful. One reason might be that the particular algorithm used only a short post-stimulus time window for estimation of the SRF, different from all other approaches; this might cause inaccuracy in the estimated SRF.
Thirdly, predictive validity could be improved by optimising high-pass filtering. A filter is commonly applied in SCR analysis to remove low-frequency drifts and condition the data. Filter choice has been based on convention or upon considerations about the response frequency spectrum (Boucsein, 2012). Here, we used an empirical approach to approximate the predictive validity of a linear method with different filter settings. We found an optimal filtering between 0.035 and 0.06 Hz with no decisive impact of filter frequency on predictive validity over this frequency range. As an intermediate choice, we suggest a 0.05 Hz unidirectional 1st order Butterworth filter as appropriate, in the context of general linear convolution modelling for eSCR.
We note that for a contrast between positive-arousing and neutral images, some aspects of the model comparison were less pronounced than for the contrasts between negative-arousing and neutral images. In particular, the precise choise of optimal filter frequency was less clear for the latter contrast. Furthermore, a non-linear model-based method -that was not informed about the experimental design -had as high a predictive validity as the best linear model. This unexpected finding differed from findings from a contrast of negative-arousing and neutral images. A possible reason is that all other methods assume a monophasic responseif responses are biphasic, these methods will be suboptimal. However, the DCM for event-related responses can easily be extended to biphasic responses as demonstrated in a previous report (Bach, Daunizeau, et al., 2010). We hope to look at this possibility in a future experiment that is optimised to assess response shape. However, no method performed better for the contrast of positivearousing and neutral images than the best linear model, and the overall pattern of results was the same for all contrasts. The fact that results are slightly less clear for this particular contrast might also reflect that the assumption of a difference in sympathetic arousal between positive and neutral images is not as unambiguously supported as for negative and neutral images (where all methods captured this difference, to a variable extent): some studies using operational approaches have suggested that sympathetic arousal in response to positive imagery might be restricted to its particular contents (e.g. explicit pictures which were not used in the present study) (Bernat, Patrick, Benning, & Tellegen, 2006).
Note we did not test the ability of the various methods to detect responses within one condition, as in a one-sample t-test on estimated SNA amplitudes for one condition. This is because such a comparison is biased and will favour methods that have skewed null distributions such as peak-scoring, model-based methods using reconstructed response functions, or non-linear models where the neural response is specified to be positive.
Finally, linear models are better able to recover SN amplitudes on a trial-by-trial basis than a non-linear DCM, although a linear, but iterative DCM for evoked responses was even better, as indicated by higher proportion of explained variance in the amplitude estimates. Note that generally the proportion of explained variance in singletrial SNA estimates is relatively low (around 0.15), reflecting large within-condition variability.
In summary, we have further validated a previously proposed linear convolution model approach to estimating SNA evoked by short arousing events. We were able to improve the scheme by augmenting the peripheral response function basis set and the filtering settings. The improved algorithm will be available as standard method in future releases of the software SCRalyze, available at http://scralyze.sourceforge.net.