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; Andé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 ψ(l)λ 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
ψ (l) λ =
1 λ ψ (l) 0 ( t λ ) ⇐⇒ ψ̂(l)λ = ψ̂ (l) 0 (λω) (1)
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:∫
ψ (l) λ (t)dt = 0 ⇐⇒ ψ̂ (l) λ (0) = 0,∀λ ∈ Λ (l), (2)
which has the following equivalent sufficient condition on each mother wavelet∫ ψ
(l) 0 (t)dt = 0 ⇐⇒ ψ̂ (l) 0 (0) = 0. (3)
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
Λ(l) = {2i/Q (l)
|i = 0, ..., J (l) ×Q(l)}. (4) 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 lth 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 LittlewoodPaley 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:
x(t) = ∑ λ1∈Λ1 (X(0)[x] ∗ ψ(1)λ1 )(t). (6)
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∫
φ(l)(t)dt = 1 ⇐⇒ φ̂(l)(0) = 1,∀l. (7)
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
φ̂(l)(ω) = e − ω2 2σ(l)2 . (8)
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 lth 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
P[c] = |c|2, ∀c ∈ C. (10)
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
X0[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]||22, 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 ofX (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
card ( {s ∈ L2(C)|S = (k1, ..., kn)} ) ≥
card ( {s ∈ L2(C)|S = (k1, ..., kn) and V = (p1, ..., pn)} ) ,
(14)
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
supp [ψ (l) λ ] := {ω|ψ (l) λ (ω) > , ω ∈ [0, 2π]} (15)
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
ψ̂µ0,σ0(ω) = H(ω)e − (ω−µ0)
2
2σ20 , (16)
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
µ0 = π 2 (2 −1/Q + 1)
σ0 = √ 3(1− 2−1/Q). (17)
Yet, instead of using the definition of scaling as defined in section 1.2.1 we will use these two parameters as follows
ψ̂µ0,σ0(λω) = ψ̂µ0λ , σ0 λ (ω) := ψ̂λ(ω). (18)
In fact, we have the following relation between the scale and the mother hyper-parameters
ψ̂µ0,σ0(λω) =H(λω)e − (λω−µ0)
2
2σ20
=H(ω)e − λ2(ω−µ0 λ )2 2σ20
=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:
supp [ψ (l) λ ] = [ µ0 λ − σ0 λ √ −2 log( ),
µ0 λ + σ0 λ
√ −2 log( ) ] .
(19)
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:
supp [φ (l)] = [ −σ(l) √ −2 log( ), σ(l) √ −2 log( ) ] . (20)
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
F [|x|2] = F [x] ∗ F [x]∗ (21)
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
2dlog2(2M)e, (22)
where dlog2(2M)e 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.
F [|x|2] = F−1 [ F [F [x]] ⊙ F [F [x]]∗ ] , (23)
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 computeX (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
S(l)λ1,...,λl [x] = 1
N N∑ t=1 X (l)λ1,...,λl [x](t). (24)
In this case, the corresponding result in Fourier is given by S(l)[x] = F [ X (l)[x] ] (ω = 0). (25)
The invariant dispersion coefficients are extracted from the Fourier transform in a straightforward manner as shown in the Appendix which results in
||X (l)[x]− S(l)[x]||22 = ||F [ X (l)[x] ]⊙ (1−F [ φ(l) ] )||22. (26)
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
(1−F [ φ(l) ] )(ω) = δ(ω) (27)
where δ denotes the Dirac function. As a result, the dispersion coefficients can be calculated as
V(l)[x] = 2 π∑
ω=1/N
( F [ X (l)[x] ] (ω) )2 , (28)
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
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.
ACKNOWLEDGEMENT
We thank Institut Universitaire de France for the Glotin’s Chair in ’Scene Analysis’. We thank SABIOD.ORG Mission Interdisciplinaire of the CNRS and Alexis Joly for co-organisation of the LifeClef Bird Challenge. We also want to thank Mr. Romain Cosentino for his reviewing work and his help in bringing back in the physical domain some of our original sentences.
A NONLINEAR INVARIANT IN THE FOURIER DOMAIN
||X (l)[x]− S(l)[x]||22 = ∫ ( X (l)[x](t)− S(l)[x](t) ) ( X (l)[x](t)− S(l)[x](t) )∗ dt
= ∫ g(t)g(t)∗dt g(t) = X (l)[x](t)− S(l)[x](t)
= ∫ F [g](ω)F [g∗](ω)dω Plancherel Theorem
=||F [g]||22 =||F [X (l)[x]− S(l)[x]]||22 =||F [ X (l)[x] ] −F [ S(l)[x] ] ||22 Linear Operator
=||F [ X (l)[x] ] −F [ X (l)[x] ∗ φ(l) ] ||22
=||F [ X (l)[x] ] −
F [ F−1 [ F [X (l)[x]] ⊙ F [φ(l)] ]] ||22
=||F [ X (l)[x] ] −F [ X (l)[x] ]⊙ F [ φ(l) ] ||22
=||F [ X (l)[x] ]⊙ (1−F [ φ(l) ] )||22.
B ADDITIONAL MATERIAL AND BIRD SONG REPRESENTATIONS
Using these three examples, we also present in Fig. 13 the resulting features computed on the first two layers of the scattering network in order to highlight the possibly linear hyperplanes separating these 3 species in this new feature space.
I find the general direction of the work is promising but, in my opinion, the paper has three main drawback. While the motivation and overall idea seem very reasonable, the derivation is not convincing mathematically. The experiments are limited and the presentation needs significant improvement.  The writing and wording are in general poorly structured to the point that it is sometimes difficult to follow the proposed ideas. The overall organization needs improvement and the connection between sections is not properly established. The paper could be significantly improved by simply re-writing it.
I'm not fully convinced by the motivation for the proposed non-linearity (|c|^2), as described on page 5. The authors argue that  (Waldspurger, 2016) suggests that higher order nonlinearities might be  beneficial for sparsity. But unless I'm missing something, that work seems  to suggest that in the general case higher order nonlinearities can be neglected. Could you please comment on this?
On the other hand, adding a second order term to the descriptor seems
an interesting direction, as long as stability to small variations is preserved (which should be shown experimentally)
The experimental section is rather limited. The paper would be stronger with a thorough numerical evaluation. The presented results, in my opinion, do not show convincingly a clear advantage of the proposed method over a standard implementation of the scattering transform. In order to show the merits of the proposed approach, it would be really helpful to directly compare running times and compression rates.
Questions: 
- Can you show empirically that the proposed higher order nonlinearity
produces sparser representations than the complex modulus?
Other minor issues:
- The proof of Section 2.1, should be preceded by a clear statement in the form of a proposition
- "Hadamart" -> Hadamard
- "Valid set" -> Validation set
- "nonzeros coefficients" -> nonzero coefficients
- Figure 3 is difficult to understand. Please provide more details.
- Figure 5 is supposed to show a comparison to a standard implementation of the Scattering network, but it doesn't seem to be such comparison in that figure. Please explain.
- Please verify the references. The first reference states "MALLAT".