Title: LINEAR TIME COMPLEXITY DEEP FOURIER SCATTERING NETWORK AND EXTENSION TO NONLINEAR INVARIANTS

Abstract: In this paper we propose a scalable version of a state-of-the-art deterministic timeinvariant feature extraction approach based on consecutive changes of basis and nonlinearities, namely, the scattering network. The first focus of the paper is to extend the scattering network to allow the use of higher order nonlinearities as well as extracting nonlinear and Fourier based statistics leading to the required invariants of any inherently structured input. In order to reach fast convolutions and to leverage the intrinsic structure of wavelets, we derive our complete model in the Fourier domain. In addition of providing fast computations, we are now able to exploit sparse matrices due to extremely high sparsity well localized in the Fourier domain. As a result, we are able to reach a true linear time complexity with inputs in the Fourier domain allowing fast and energy efficient solutions to machine learning tasks. Validation of the features and computational results will be presented through the use of these invariant coefficients to perform classification on audio recordings of bird songs captured in multiple different soundscapes. In the end, the applicability of the presented solutions to deep artificial neural networks is discussed.

1. INTRODUCTION AND SCATTERING NETWORK:

1.1 BACKGROUND Invariants are the gems of machine learning enabling key latent space representations of given inputs. Following this analogy, precious invariants shine out by being discriminative enough to detect changes in the underlying data distribution yet with bounded variations to ensure stable and robust representations. The motivation to find informative invariants as latent representations is becoming the principal focus from deep learning to signal processing communities aiming to tackle most machine learning tasks. Undoubtedly, given infinite datasets and computational power, learning invariants will lead to the fittest descriptors. However, nowadays problems do not fit this situation forcing the use of nested architectures and supercomputers to reach, in the limit, these utopian descriptors. As an alternative, the scattering network (Mallat, 2012; Bruna & Mallat, 2013; Andén & Mallat, 2014) provides a deterministic transformation of a given input signal x through a cascade of linear and nonlinear operators which do not commute. The linear operator is able via a dimensional increase to linearize the representation in the sense that x can be expressed as a linear combination of basic yet fundamental structures. This linear transformation is in practice a wavelet transform but can be generalized to any complete or over-complete change of basis. Originally, these wavelet transforms were used with an over-complete basis derived from Morlet and Gabor wavelets (Mallat, 1989) . Recently, a discrete wavelet transform scheme (Mallat, 1999) and specifically a Haar transform (Chen et al., 2014) has been used instead to reduce the computational overload of the scattering transform. This framework, however, is not suited for general tasks due to poor frequency resolution of one wavelet per octave and the not continuously differentiable Haar wavelet (Graps, 1995) making it unsuitable for biological and natural waveforms detection. Following the change of basis, a k-Lipschitz nonlinear operator is applied which must also be contractive to enforce space contraction and thus bound the output variations (Mallat, 2016) . The nonlinearity used in the scattering network is the complex modulus which is piecewise linear. This surjection aims to map the transformation into a subspace of smaller radius {|x| | x ∈ Ω} ⊂ Ω where Ω is the space of studied signals. As a result, one can see these successions of transforms as a suite of expansion and contraction of a deeper and deeper signal representation in the hope to decode, detect and separate all the underlying energy sources. These layered representations however still contain the time dimension and are thus overall extremely sparse and not translation invariant. This time sensitivity motivates the aggregation of the time dimension. This led to the scattering coefficients per se which are computed through the application of the operator S applied on each previously computed representation and each frequency band. It is defined as an order one statistic, the arithmetic mean, over the time dimension leading to a time-invariant central tendency description of the layered representation of x. The resulting scattering coefficients, when used as features in machine learning tasks, led to state-of-the-art results in music genre classification (Chen & Ramadge, 2013) , texture classification (Sifre & Mallat, 2013; Bruna & Mallat, 2013) and object classification (Oyallon & Mallat, 2015) . In this paper, we present a modification of the standard scattering network by replacing the complex modulus nonlinearity with a quadratic nonlinearity in order to increase the SNR while allowing to compute the complete scattering network without leaving the Fourier domain, which was before necessary after each level of decomposition.

1.2. SCATTERING NETWORK:

We now present formally all the steps involved in the scattering network in order to clarify notations and concepts while making it as explicit as possible to easily develop our extensions. For readers more familiar with neural networks, it is first important to note that this framework can be seen as a restricted case of a Convolutional Neural Network (LeCun & Bengio, 1995) where the filters are fixed wavelets and the nonlinearity is the complex modulus as well as some topological differences such as depicted in Fig. 1 . The scattering coefficients are then computed through time averaging of each representation.

1.2.1. HIERARCHICAL REPRESENTATION:

By definition a scattering network can have any fixed number of layers L defined a priori. These layers are ordered in a hierarchical fashion so that the output of layer l is the input of layer l + 1. In the following, many presented properties and definitions hold for all l ∈ {1, ..., L}. Each layer l uses a specific filter-bank made of band-pass filters ψ 

EQUATION

The collection of scaling factors λ for layer l is denoted by Λ (l) . The only admissibility condition that each filter must satisfy is to have zero mean:

EQUATION

which has the following equivalent sufficient condition on each mother wavelet

EQUATION

The finite set of continuous scale factors needed to derive the filer-bank is given as a geometric progression governed by two hyper-parameters, the number of wavelets per octave Q and the number of octave to decompose J. The Q parameter, also called quality criteria, defines the frequency resolution, the greater it is the finer the resolution is but the more redundant will be the representation.

The J parameter defines the number of octave to decompose. Since these parameters can be layer specific we now denote them as Q (l) and J (l) . We thus have

EQUATION

When the L filter-banks are generated, it is possible to compute the L representations by iteratively applying the filter-banks and the nonlinearity. As a result, the l th representation indexed by the time dimension t with the l first scales as hyper-parameters is given by:

X (0) [x](t) := x(t) X (1) λ1 [x](t) := |(X (0) [x] * ψ (1) λ1 )(t)|, ∀λ 1 ∈ Λ 1 X (l) λ1,...,λ l [x](t) = |(X (l-1) λ1,...,λ l-1 [x] * ψ (l) λ l )(t)|, ∀λ 1 ∈ Λ (1) , ..., λ l ∈ Λ (l) (5)

One can notice that X

(1) λ1 [x] coefficients form the well known wavelet transform or scalogram. We now denote by X (l) [x] the complete time-frequency representation for layer l if the scales are not necessary in the context.

Thinking of deeper layers as representations of more abstract concepts is inappropriate and thus not analogous to deep neural networks representations simply because deeper layer filters are not linear combination of the first layer filters. Since the used filters are renormalized to satisfy the Littlewood-Paley condition, the energy contained in each layer decays exponentially (Waldspurger, 2016) . As a result, deeper layer will contain less and less energy until all events have been captured and all the next layers are zeros. With this renormalization, inverting one change of basis is instantaneous, simply add up together all the coefficients obtained from the filters application:

EQUATION

This last property highlights one motivation of dimensional increase, event or structure separation.

It is now possible to rewrite the treated signal as a combination of fundamental structures, namely, the responses of the signal with the filters linearizing the events w.r.t the new λ dimension.

1.2.2. SCATTERING COEFFICIENTS:

For each of the L representations, one can extract the scattering coefficients by applying a scaling function φ (l) on X (l) [x] leading to the time invariant representation S (l) [x]. The scaling function acts as a smoothing function on a given time support and satisfies

EQUATION

The scaling function is usually a Gaussian filter with layer dependent standard deviations 1/σ (l) in the time domain and σ (l) in the frequency domain defined as

EQUATION

)

The greater the standard deviation is in the physical domain the more time invariant are the scattering coefficients. Ultimately, we reach global time-invariance and the scattering operator S[x] reduces to an arithmetic mean over the input support. Since only the standard deviation of the scaling function is layer dependent, we denote by φ (l) the l th scaling function generated using σ (l) . We can thus define the scattering operators as

S (0) [x](t) := (X (0) [x] * φ (0) )(t), S (1) λ1 [x](t) := (X (1) λ1 [x] * φ (1) )(t), ∀λ 1 ∈ Λ (1) S (l) λ1,...,λ l [x](t) := (X (l) λ1,...,λ l [x] * φ (l) )(t) ∀λ 1 ∈ Λ (1) , ..., λ l ∈ Λ (l) (9)

For machine learning tasks, using global time-invariance yield robust yet biased time invariant descriptors due to too much many-to-one possible mappings. As a result, a local or windowed scattering transform has been used (Bruna & Mallat, 2013) leading to only local time-invariance through a smaller time support for the scaling function. This tweak is possible in computer vision tasks where each input is of the same size and the role of φ is to bring local diffeomorphism invariance which has been shown to smooth the underlying manifold (Zoran & Weiss, 2011; 2012) . However, for audio tasks and more general problems, this constant input size is rare forcing the use of φ for completely aggregating the time dimension.

2. EXTENSIONS:



2.1. HIGHER ORDER NONLINEARITY:

The usual nonlinearity applied in a scattering network is the complex modulus. This nonlinearity is not everywhere differentiable but is contractive leading to an exponential decay in the energy distribution over the layers. However, as pointed out in (Waldspurger, 2016) , higher order nonlinearity might be beneficial sparsity-wise and to increase the SNR. As a result, we chose to use a continuously differentiable second order nonlinearity which will have the beneficial property of adapting its contractive property for irrelevant inputs while maintaining bounded variations. This nonlinearity is defined as

EQUATION

)

Proof. We now prove the adaptive k-Lipschitz property of our nonlinearity

||P[a] -P[b]|| =|||a| 2 -|b| 2 || =||(|a| -|b|)(|a| + |b|)|| =(|a| + |b|)|||a| -|b||| ≤(|a| + |b|)||a -b|| =K(a, b)||a -b||

Since the input signal is renormalized so that ||x|| 1 = 1, we have that |a| + |b| ∈ [0, 1[. As a result, given the inputs constraints, P is a contractive operator with bounded variations. Yet, the degree of contraction will vary given the input amplitudes leading to a better SNR. This means that high amplitudes resulting from close match between the filter and the signal will be efficiently represented whereas small amplitude coefficients resulting from noise filtering and mismatches between the filter and the signal will be highly contracted. Since in practice wavelet filters catch the relevant events, this property allows high quality representations. This change will not only increase the relative sparsity of the representations but also allow us to perform some major computational tricks as describe in the following section. As a result we define the new representations as

X 0 [x](t) := x(t) X (1) λ1 [x](t) := P (X (0) [x] * ψ (1) λ1 )(t) , ∀λ 1 ∈ Λ (1) X (l) λ1,...,λ l [x](t) = P (X (l-1) λ1,...,λ l [x] * ψ (l) λ l )(t) , ∀λ 1 ∈ Λ (1) , ..., λ l ∈ Λ (l) (11)

as well as the new scattering coefficients as

S (0) [x](t) := (X (0) [x] * φ (0) )(t), S (1) λ1 [x](t) := (X (1) λ1 [x] * φ (1) )(t), ∀λ 1 ∈ Λ (1) S (l) λ1,...,λ l [x](t) := (X (l) λ1,...,λ l [x] * φ (l) )(t), ∀λ 1 ∈ Λ (1) , ..., λ l ∈ Λ (l) (12)

2.2. INVARIANT DISPERSION COEFFICIENTS:

The scattering coefficients used to characterize the signal of interest are known to be efficient for stationary inputs but not descriptive enough otherwise. We thus propose to generate complementary invariant coefficients based on a dispersion measure, the variance. As a result, these complementary coefficients derived from the second order moment will help to characterize the input leading to more discriminative features while maintaining global time invariance. We now define these invariant dispersion coefficients as

V (0) [x] := ||X (0) [x] -S (0) [x]|| 2 2 , V (1) λ1 [x] := ||X (1) λ1 [x] -S (1) λ1 [x]|| 2 2 , ∀λ 1 ∈ Λ (1) V (l) λ1,...,λ l [x] := ||X (l) λ1,...,λ l [x] -S (l) λ1,...,λ l [x]|| 2 2 . ∀λ 1 ∈ Λ (1) , ..., λ l ∈ Λ (l) . (13)

The resulting V (l) [x] coefficients are thus globally time invariant whatever scaling function was used to compute S (l) [x]. In fact, these invariant dispersion coefficients represent the variance between X (l) [x] and S (l) [x] representations whether S (l) [x] was globally time invariant or a smoothed version of X (l) [x]. In order for V (l) [x] to be invariant to random permutations as well, S (l) [x] should be globally translation invariant and thus also globally invariant to random permutations. In addition, regarding the discriminative ability gained through the use of these second order statistic, we have that

EQUATION

where S represents a realization of the scattering coefficients for all layers, all frequency bands, and V a realization of the dispersion coefficients again for all layers and all frequency bands. From this, it follows that the set of invariant coefficients

(S[x], V[x]

) is more discriminative leading to more precise data description than when using (S[x]) only. The development of these presented invariant dispersion coefficients opened the door to the development of uncountably many new invariant coefficients. We now present the elaboration of the scheme in the Fourier domain and the computational tricks involved in order to reach linear complexity.

3. FAST IMPLEMENTATION, FOURIER DOMAIN AND LINEAR COMPLEXITY:



3.1. INTRODUCTION AND SPARSE STORAGE:

One of the great benefits of the wavelet transform is the induced sparsity in the representation for certain class of signals (Elad & Aharon, 2006; Starck et al., 2010) which is seen as a quality criteria of the representation (Coifman et al., 1992) . In addition of providing sparse representations, wavelets are localized filters in time and frequency domain. However, the idea to exploit this known sparsity in order to reduce the computational time of performing a transformation has not been leveraged yet. When dealing with the standard time domain, one can not know a priori where the filter with match or not the signal and thus where are the nonzeros coefficients leading to no way to have computational gains. On the other hand, applying the filter in the Fourier domain reduces to an Hadamard product and thus the resulting support is deduced from the filter support which is known to be localized. As a result, it is now possible to know a priori most of the zero coefficients positions since in the Fourier domain the filter is well localized. Also, the Fourier domain, the wavelet support is convex and compact but most importantly it can be computed a priori given the scale parameter and the mother wavelet. This motivates our choice to perform our framework including the wavelet transform, the nonlinearity P and the invariant features extraction in the Fourier domain leading to linear complexity. Furthermore, using the Fourier domain allows us to efficiently leverage sparse matrices leading to efficient storage and memory management on energy efficient platforms such as presented in (Esser et al., 2015) . We will first present the computation of the filters in the Fourier domain as well as their convex compact support derivation. From that we present the sparse application of the filters and how to compute the nonlinearity in Fourier. Finally, we will see that extracting the invariant features can be done efficiently leading to our main result which is a linear complexity overall framework. Concerning the Fourier transform, the Danielson-Lanczos lemma (Flannery et al., 1992) will be used in order to provide a true O(N log(N )) complexity for an input of size N which is a power of 2. As we will see, this requirement will always be fulfilled without any additional cost.

3.2. SPARSE FOURIER FILTERS:

One particularity of the continuous wavelets such as DoG, Morlet wavelets reside in their localized compact support in the Fourier domain. In our description the used wavelet will be a Morlet wavelet but this analysis can be extended to any continuous wavelet with analytical form. We define the support of the filter ψ (l)

λ given the threshold as

EQUATION

As presented in (Balestriero et al., 2015) the scales define entirely the support of each wavelet. In order to develop synergistic computational tricks, we derive our framework in the Fourier domain.

Let define the Morlet wavelet as

EQUATION

)

where the parameters µ 0 and σ 0 represent respectively the center frequency and bandwidth of the mother wavelet and H is the step-wise function. The ratio between these two quantities will remain the same among all the filters, in fact, wavelets have a constant ratio of bandwidth to center frequency. These two mother hyper-parameters are taken as

EQUATION

Yet, instead of using the definition of scaling as defined in section 1.2.1 we will use these two parameters as follows

EQUATION

)

Figure 2 : Filter-Bank in Fourier with 4 wavelets per octave on 4 octaves. The sparsity is of 0.94% with = 0.01 ) needed to compute the filter-bank given the signal size with standard parameters J = 5, Q = 16 growing linearly with the input size on 1 CPU @1.4GHz.

Signal Size 2 16 2 18 2 20 2 22 Time in sec. 0.006 0.018 0.066 0.262

In fact, we have the following relation between the scale and the mother hyper-parameters ψµ0,σ0 (λω) =H(λω)e

- (λω-µ 0 ) 2 2σ 2 0 =H(ω)e - λ 2 (ω- µ 0 λ ) 2 2σ 2 0 =H(ω)e - (ω- µ 0 λ ) 2 2( σ 0 λ ) 2 . = ψ µ 0 λ , σ 0 λ (ω)

Given this, we can compute explicitly the support of every filter ψ (l)

λ . As one can notice these filters have a convex compact support around their center frequencies:

EQUATION

In Fig. 2 one can see a filter-bank example where all the filters are presented in one plot demonstrating the important sparsity inherent to wavelets. Varying the parameter affects directly the number of nonzero coefficients and we thus also present in Table 1 the exact sparsity with different input sizes and parameter. Since the support is known a priori, it is straightforward to optimize their computation and allocation through sparse vectors leading to fast filter-bank generation as presented in Table 2 with large input sizes. In fact, the input length defines the size of the generated filter since we now perform the convolution in the Fourier domain.

Concerning the φ (l) filter, given its bandwidth σ (l) , its support is given by:

EQUATION

Some examples are shown in Fig. 8 where different σ (l) are selected representing different time supports. For each of the filters and each of the layers, the application is done through element-wise multiplication in the Fourier domain as explained in the next section where the nonlinearity will be defined.

3.3. NONLINEARITY AND FILTERING IN FOURIER:

The nonlinearity used in this framework defined in section 2.1 is efficiently done in the Fourier domain through the following property

EQUATION

If done directly, this operation would be slower in the Fourier domain since we jump from a linear complexity to a quadratic complexity. However, one should notice from section 3.2 that we are dealing with F[x] which are extremely sparse but most importantly with convex compact support of size M << N . Exploiting this sparsity could lead to a faster convolution which would still be of quadratic complexity w.r.t the support size. However, using the convolution theorem it is possible to perform this convolution in M log(M ) complexity by applying again a Fourier transform now only on the convex compact support of F[x]. In order to have proper boundary condition and not the periodic assumption of the Fourier transform we use a zero-padded version of size2M instead of M leading to exact computation of the convolution. In addition, the support size of 2M is the minimum required size. For a fast Fourier transform algorithm, this has to be a power of 2. As a result, in practice, the zero padding is done to reach the size which is the smallest power of 2 greater than 2M defined as

EQUATION

where log 2 (2M ) denotes the smallest of the greater integers. As a result in the Fourier domain we will apply another Fourier transform in order to compute this auto-correlation which will correspond to the desired nonlinearity in the time domain.

EQUATION

)

where is the Hadamard product, F is the Discrete Fourier Transform and F -1 its inverse operator. In addition of the second Fourier transform being applied on a really small support, it is also important to note that after application of the nonlinearity the output is conjugate symmetric in the Fourier domain but since the filter-banks are always applied on [0, π] we can store only this part for further computation and re-generate the conjugate symmetric part when applying φ (l) . We present this operation in the Fourier domain in Fig. 7 .

In order to highlight the high sparsity encountered in the Fourier domain when dealing with this filter application, we present in Fig. 3 an example where the nonzero elements are shown. This corresponds to the first representation namely X (1) [x]. For the second representation, the input support will not be over the whole [0, 2π[ domain but around 0 and thus implies increased sparsity as demonstrated in Fig. 9 . Given these two descriptions, one is able to compute X (l) [x] for any l. We thus now present how to compute the scattering and dispersion coefficients given this representations in the Fourier domain.

3.4. SCATTERING COEFFICIENTS EXTRACTION:

The scattering coefficients S (l) [x] result from the application of a Gaussian filter parameterized by its standard deviation. In the general case where global time invariance is required, this standard deviation is taken to be infinite in the time domain resulting in

EQUATION

)

Figure 3 : Nonzero elements are shown in black, after application of the filters and the nonlinearity operator P. The x-axis corresponds to the frequency ω and the y-axis to the filter index. This follows the geometric progression of the wavelets since each line corresponds to ψ(i) x.

In this case, the corresponding result in Fourier is given by

EQUATION

)

The invariant dispersion coefficients are extracted from the Fourier transform in a straightforward manner as shown in the Appendix which results in

EQUATION

Thus (1-F φ (l) ) acts as a mask to reduce the norm computation by the amount of energy captured through the scaling function application. For the case where we have global time invariance or infinite standard deviation, this mask reduces to

EQUATION

where δ denotes the Dirac function. As a result, the dispersion coefficients can be calculated as

EQUATION

which is the L2 norm computed without taking into account the coefficient at ω = 0 exploiting the conjugate symmetry structure for the real input signal x. Conceptually, the V coefficients capture the remaining energy and thus ensures that for any depth of the scattering network, all the energy of the input is contained inside the computed invariants. In fact, one can see that

V (l) [x] = ∞ i=l+1 S (i) [x].

3.5. SCALABILITY:

We now present some results and figures in order to highlight the high scalability of our framework with respect to the input size. Note that the current implementation is done in Python. Implementing this toolbox in C is a future work which will lead to even better results than the ones displayed below which are nevertheless already astonishing. First of all, one can see in Fig. 4 that the number of nonzero coefficients increase linearly with the input size. This result is important in nowadays paradigm where technologies allow extreme frequency sampling and thus input signals with high dimensions yet we aim to save as much memory and storage as possible. If we put this nonzero coefficients in perspective with the total possible number of coefficients we obtain our sparsity coefficient which grows logarithmically with the input size as shown in Fig. 4 . This result shows the advantage of using sparse matrices which increases as the input size increases. The sparsity is thus in our case a justification to exploit the Fourier domain.

Finally, in Fig. 5 are presented some computational time for different input signals. We can see in this figure the high efficiency of our approach put in perspective of an existing C implementation of the scattering network. In fact, in this latter, one had to perform multiple inverse Fourier transforms in order to apply the nonlinearity and in order to compute the second layer for example apply again a Fourier transform and this for all the frequency bands. As a result the previously fastest known algorithm was of asymptotic complexity O(N log(N )) even with a Fourier input. In addition, it did not leverage the sparsity of the representation leading to poor memory management and storing. Not however that for all the existing implementations, the complexity is linear with respect to the J and Q parameters. Finally, with our framework, one can directly store the sparse matrices of the representations leading to space saving on the hard drive in addition of the actual Random Access Memory (RAM) saving during the computation.

4. VALIDATION ON BIRD CHALLENGE:



4.1. DATASET PRESENTATION:

We now validate the use of the scattering coefficients as features for signal characterization through a supervised problem of audio recordings classification. The bird song classification challenge is made of 50 classes and correspond to a small version of the BirdCLEF Challenge (Joly et al., 2015) . The recordings come from the Xeno-Canto database and mainly focus on species form South America.

For our task, the dataset used for training is composed of 924 files sampled at 44100 Hz. The validation set used to present our classification accuracy contains about 400 files. Computing the S[x] and V[x] features on the training and validation set takes between 2 to 3 hours depending on the set of parameters used with a 2-layer architecture on 1 CPU. The files add up to a disk usage of 4.2Go, the computed set of features however represent 450Mo. As a result, we are able to encode and extract important characteristics of the signals while reducing the amount of redundant information. We present in Fig. 101112 examples of the dataset with the waveform as well as the representation X

(1)

λ1 [x].

The aim is to first demonstrate the sparsity or high SNR in the physical domain involved by using a second order nonlinearity. In addition, one can see the different frequency modulated responses that could characterize a bird song. Overall, there are some fundamental difficulties in this task. The first challenge is to deal with translation invariance. In fact, the bird songs can be captured anywhere inside each files which themselves are of many different durations, from seconds to minutes. The second difficulty resides in characterizing well enough the time-frequency patterns of each specie without being altered by the ambient noise or possible presence of other species including human voices. Finally, difficulties also arise from the machine learning point of view with large class imbalance in the dataset.

4.2. RESULTS:

Figure 6 : Classification MAP given two set of parameters for the second scattering network layer. Are also presented the results when using only the scattering coefficients S[x], the dispersion coefficients V[x] and a concatenation of both, best is 52%.

We now present the classification results obtained via our developed framework. First of all, no additional features have been engineered and thus only the features developed in the paper are used. For the classification part, we decided to use a fast algorithm in accordance with the whole scheme developed earlier and thus used random forests (Liaw & Wiener, 2002) . In short, random forests are using bagging of decision trees (Breiman, 1996) and thus are able to aggregate multiple weak classifiers to end up with efficient class boundaries. One of its drawback resides in the fact that it can only create decision rule on each feature dimension without combining them as could do a logistic regression for example. In addition, we used a weighted loss function in order to deal with the imbalanced dataset (Van Hulse et al., 2007) . Finally, no additional pre-processing/denoising has been used and no feature extraction/selection technique has been piped in. Yet, with this basic approach, we were able to reach an accuracy of 47.6% and a Mean Average Precision (MAP) of 52.4%. The state-of-the-art technique for this problem reached a MAP of 53% (Cha, 2016). We present in Fig. 6 some accuracy results where two sets of parameters have been used for the second layer of the scattering network. In addition, we show the classification results when using each features independently and combination of the two in order to highlight their complementarity. Given the deterministic transformation used and the lack of cross-validation and fine tuning, we think of these results as promising overall while being state-of-the-art if considering solutions where no learning was involved outside of the classifier. For example, one extension on the classifier could be to use boosting algorithms (Schapire et al., 1998) or neural networks. Concerning the representation, performing cross-validation on the parameters could lead to great improvements as finally a third scattering layered could also be considered.

5. CONCLUSION:

We presented an extension of the scattering network in order to provide more discriminative time invariant features which are complementary to each others. The derivation of a second order invariant operator as well as the use of a second order nonlinearity in the layered representation computation led to efficient characterization of audio signals opening the door of many more possible time invariant features derivation. The whole framework has been derived in the Fourier domain in order to reach linear complexity in the input size as it is now possible to compute all the layer representations and feature without leaving the Fourier domain. Sparse storage is also a milestone of our algorithm leading to not only efficient computation but smart memory management allowing this framework to be applied online on energy efficient laptops and chips. In addition, only simple arithmetic operations are used and parallel implementation can be done easily as well as GPU portage. This framework can be applied without any assumption on the input signal and thus aims to be as general as possible as a unsupervised invariant feature extraction. Finally, we hope to bring the consideration of sparse filters and Fourier based computation for deep convolutional networks. In fact, as the datasets get larger and larger, the complexity of the networks increase and convolutions might not be efficiently computed in the physical domain anymore. Since the convergence of the filter ensure their sparsity and smoothness, this consideration might help to bring deep learning to the family of scalable algorithms with the development of Fourier networks as a whole.

Figure 7 : Demonstration of the P computation in the Fourier domain. Top: input after application of a specific filter. Middle: Extracted window with nonzero elements and optimal padding greater than the minimum size up to the next power of 2. Bottom: Result of the convolution done through another Fourier transform and the convolution theorem, the kept coefficients are from 0 to M since they are followed by zeros and the complex conjugate of these coefficients leading to optimal results. 

:



Figures and Tables:

derived by scaling the mother wavelet of layer l denoted as ψ (l) 0 in the time domain and ψ(l) 0 in the Fourier domain. The dilation factors are denoted by the subscript λ leading to

Figure 1: Architecture difference between the CNN (A) and the scattering network (B) without depiction of the computation of the scattering coefficients which are obtained after averaging over each obtained representation.

Figure 4: Left: Figure showing the increase in nonzero coefficients is linear with respect to the input size. Right: Figure showing the increase in sparsity in our representations for the two layers. The increase is logarithmic with respect to the input size.

Figure 5: Needed computation time to perform the transform for the first X (1) [x] and second X (2) [x] layers. The need computation time is more than 16 times smaller than the known Scatnet toolbox implemented in C (github.com/RandallBalestriero/CIGAL_GUI) of O(N log(N )) time complexity even when dealing with Fourier domain inputs. This shows the advantage of our approach which here is implemented in Python only.

Figure 8: Some possible φ filters in the Fourier domain corresponding to Gaussian filtering with bandwidth in the physical domain inversely proportional to the σ in the Fourier domain without renormalization.

Figure 9: nonzero elements present after application of the filters and the nonlinearity operator P on this representation for a sparse input.

Figure 10: Example 1: transform X (1) λ1 [x]. In this case, clear frequency modulations appear for only one source and high SNR. The noise is contracted to 0 through the nonlinearity.

Figure 11: Example 2: transform X (1) λ1 [x]. In this case, the source presents multiple kinds of chirps, and frequency modulated patterns. Some harmonics are detected yet it is clear that aggregation of the time dimension with this representation only will aggregate the different patterns leading to poor signal characterization leading to the need of X (2) λ1,λ2 [x].

Figure 12: Example 3: transform X (1) λ1 [x]. In this case only transients are present and almost not frequency modulation appear on the features. This kind of signal will be well captured with one layer only.

Figure 13: We present here the features extracted form the 3 examples presented in Fig. 10,11,12.The left part contains the scattering coefficients encoding the arithmetic mean whereas the right part concerns the dispersion coefficients. On the top part the features are extracted from the first layer and on the bottom are the features extracted form the second layer. It is clear that for these signals, the features of the first layer are enough to discriminate them. Notice that through global time invariance, one ends up with features vectors of exact same dimension for each signal and that they would remain the same if the input signal was translated.



Table content: <table><tr><td>= 0.0001</td><td>= 0.0000001</td></tr><tr><td>N = 524288 (2 19 ) 98.39883%</td><td>97.89803%</td></tr><tr><td>N = 1048576 (2 20 ) 98.39893%</td><td>97.89812%</td></tr><tr><td>N = 2097152 (2 21 ) 98.39898%</td><td>97.89816%</td></tr></table>

Time (in sec.

Table content: <table/>

Peer review (iclr2017):
Overview: This work seems very promising, but I believe it should be compared with more baselines, and more precisely described and explained, from a signal processing point of view. Pros: New descriptor Fast implementation Cons: a) Lack of rigor b) Too long accordingly to the content c) The computational gain of the algorithm is not clear d) The work is not compared with its most obvious baseline: a scattering transform I will detail each cons. a) Section 1: The author motivates the use of scattering transform because it defines a contraction of the space that relies on geometric features. * The nonlinearity used in the scattering network is the complex modulus which is piecewise linear.* A real modulus is piecewise linear. A complex modulus has a shape of bell when interpreting C as R^2. Could you clarify? Omega is not introduced. Could you give a precise reference (page+paper) of this claim: “Higher order nonlinearity refers to |x|^2 instead of |x| as it is usually done in the scattering network.” ? Section 2: The motivation of the non-linearity is not clear. First, this non-linearity might potentially increase a lot the variance of your architecture since it depends on higher moments(up to 4). I think a fair analysis would be to compute numerically the normalized variance (e.g. divided by the averaged l^2 norm), as a sanity check. Besides, one should prove that the energy is decreasing. It is not possible to argue that this architecture is similar to a scattering transform which has precise mathematical foundations and those results are required, since the setting is different. Permutation is not a relevant variability. The notion of sparsity during the whole paper sometimes refers to the number of 0 value, either the l^1 norm. Mathematically, a small value, even 10^-1000 is still a non 0 value. Did you compute the graph of the figure 4 on the bird dataset? You might use a ratio instead for clarity. The wavelet that is defined is not a morlet wavelet ( https://en.wikipedia.org/wiki/Morlet_wavelet ). It is close to be a gabor wavelet, and actually it has not a 0 averaging. The measure of the *sparsity of the filters* is extremely unclear, is it the ratio between the support of the filter and its size? A good criteria might be for instance to understand the amount of the energy that has been neglected. Besides, a filter with compact support has a bad localisation property. However, this topic is not reached in the paper. For instance, Cauchy wavelets are not used in many applications.(However in mathematical proofs, they often are) In Subsection 3.4 you write that V^(l)(x)=Sum_{j>l} S^(j)(x), but also that you do compute only the 2 first order coefficients because they can be neglected. Besides, you specifically write that adding the variance coefficients improve the representation, whereas they can be obtained as linear combination of S. You claim you apply only one FFT, whereas you apply several FFTs. b) From *One of the great...* to *iof the input.* section 3.1, the text is not clear. The motivation is that a convolution in space is slower that performing a convolution in the Fourier domain. This whole paragraph can be summarised in few sentences. The section 3.4 is long and corresponds to implementation details. Maybe it could be removed. c) The table 2 seems to indicate that the generation of the filters is one of the bottleneck of your software. Is this really true? One of the main claim of the paper is that sparse filters and staying in Fourier domain speed up the computations. Let us compare the computation of the first order scattering coefficients at scale j with this setting. One has to compare the complexity to compute sum |x*psi_j| and sum |x*psi_j|^2 A downsampling is always performed with wavelets, yet it bears approximation. In a Fourier resolution implementation, one adjust the degree of approximation and speed of computation with the over sampling parameters. Assume a FFT of size N costs C*N*log N, then, Computing hat x costs in both case C*N*log N. ST: Then, the signal is multiplied with the fourier transform of the filter, which has a cost of N. In a fourier multi resolution implementation, one periodises the signal by a factor j, such that its size is N/2^j. Then, the FFT has cost C*N/2^j*log(N/2^j), and modulus has cost say N. Then, one applied the averaging that has complexity N/2^j. Here: The signal is multiplied with the fourier transport of the filter that has a support of N/2^j. Then, you convolve it with itself, that has thanks to the padding and the FFT a cost of C*N/2^(j-1)*log(N/2^(j-1))+N/2^(j-1). And you take the 0 frequency. I might be wrong, since this it not my work to do those calculus, but if you claim that your implementation is theoretically faster, you need to prove it, since I do not know any papers where scattering transform claims to be fastly implemented. Here, one sees that the difference is not that significant. Please correct me if I did a mistake. d) It is essential to compare your work with the representation of a scattering transform. First, in term of speed of computation, with a fair implementation (e.g. not MATLAB) and secondly in term of accuracy on the dataset you did use: it is a natural baseline.