Understanding deep convolutional networks

Deep convolutional networks provide state-of-the-art classifications and regressions results over many high-dimensional problems. We review their architecture, which scatters data with a cascade of linear filter weights and nonlinearities. A mathematical framework is introduced to analyse their properties. Computations of invariants involve multiscale contractions with wavelets, the linearization of hierarchical symmetries and sparse separations. Applications are discussed.


§1 Introduction
Supervised learning is a high-dimensional interpolation problem.We approximate a function f (x) from q training samples {x i , f (x i )} i≤q , where x is a data vector of very high dimension d.This dimension is often larger than 10 6 , for images or other large size signals.Deep convolutional neural networks have recently obtained remarkable experimental results [21].They give state of the art performances for image classification with thousands of complex classes [19], speech recognition [17], bio-medical applications [22], natural language understanding [30], and in many other domains.They are also studied as neuro-physiological models of vision [4].
Multilayer neural networks are computational learning architectures which propagate the input data across a sequence of linear operators and simple non-linearities.The properties of shallow networks, with one hidden layer, are well understood as decompositions in families of ridge functions [10].However, these approaches do not extend to networks with more layers.Deep convolutional neural networks, introduced by Le Cun [20], are implemented with linear convolutions followed by non-linearities, over typically more than 5 layers.These complex programmable machines, defined by potentially billions of filter weights, bring us to a different mathematical world.
Many researchers have pointed out that deep convolution networks are computing progressively more powerful invariants as depth increases [4,21], but relations with networks weights and non-linearities are complex.This paper aims at clarifying important principles which govern the properties of such networks, but their architecture and weights may differ with applications.We show that computations of invariants involve multiscale contractions, the linearization of hierarchical symmetries, and sparse separations.This conceptual basis is only a first step towards a full mathematical understanding of convolutional network properties.
In high dimension, x has a considerable number of parameters, which is a dimensionality curse.Sampling uniformly a volume of dimension d requires a number of samples which grows exponentially with d.In most applications, the number q of training samples rather grows linearly with d.It is possible to approximate f (x) with so few samples, only if f has some strong regularity properties allowing to ultimately reduce the dimension of the estimation.Any learning algorithm, including deep convolutional networks, thus relies on an underlying assumption of regularity.Specifying the nature of this regularity is one of the core mathematical problem.
One can try to circumvent the curse of dimensionality by reducing the variability or the dimension of x, without sacrificing the ability to approximate f (x).This is done by defining a new variable Φ(x) where Φ is a contractive operator which reduces the range of variations of x, while still separating different values of f : Φ(x) = Φ(x ) if f (x) = f (x ).This separation-contraction trade-off needs to be adjusted to the properties of f .Linearization is a strategy used in machine learning to reduce the dimension with a linear projector.A low-dimensional linear projection of x can separate the values of f if this function remains constant in the direction of a high-dimensional linear space.This is rarely the case, but one can try to find Φ(x) which linearizes high-dimensional domains where f (x) remains constant.The dimension is then reduced by applying a low-dimensional linear projector on Φ(x).Finding such a Φ is the dream of kernel learning algorithms, explained in Section 2.
Deep neural networks are more conservative.They progressively contract the space and linearize transformations along which f remains nearly constant, to preserve separation.Such directions are defined by linear operators which belong to groups of local symmetries, introduced in Section 3. To understand the difficulty to linearize the action of high-dimensional groups of operators, we begin with the groups of translations and diffeomorphisms, which deform signals.They capture essential mathematical properties that are extended to general deep network symmetries, in Section 7.
To linearize diffeomorphisms and preserve separability, Section 4 shows that we must separate the variations of x at different scales, with a wavelet transform.This is implemented with multiscale filter convolutions, which are building blocks of deep convolution filtering.General deep network architectures are introduced in Section 5.They iterate on linear operators which filter and linearly combine different channels in each network layer, followed by contractive non-linearities.
To understand how non-linear contractions interact with linear operators, Section 6 begins with simpler networks which do not recombine channels in each layer.It defines a non-linear scattering transform, introduced in [24], where wavelets have a separation and linearization role.The resulting contraction, linearization and separability properties are reviewed.We shall see that sparsity is important for separation.Section 7 extends these ideas to a more general class of deep convolutional networks.Channel combinations provide the flexibility needed to extend translations to larger groups of local symmetries adapted to f .The network is structured by factorizing groups of symmetries, in which case all linear operators are generalized convolutions.Computations are ultimately performed with filter weights, which are learned.Their relation with groups of symmetries is explained.A major issue is to preserve a separation margin across classification frontiers.Deep convolutional networks have the ability to do so, by separating network fibers which are progressively more invariant and specialized.This can give rise to invariant grandmother type neurons observed in deep networks [1].The paper studies architectures as opposed to computational learning of network weights, which is an outstanding optimization issue [21].
§2 Linearization, Projection and Separability Supervised learning computes an approximation f (x) of a function f (x) from q training samples {x i , f (x i )} i≤q , for x = (x(1), ..., x(d)) ∈ Ω.The domain Ω is a high dimensional open subset of R d , not a low-dimensional manifold.In a regression problem, f (x) takes its values in R, whereas in classification its values are class indices.
Separation Ideally, we would like to reduce the dimension of x by computing a low dimensional vector Φ(x) such that one can write We then say that Φ separates f .For regression problems, to guarantee that f 0 is regular, we further impose that the separation is Lipschitz: It implies that 2 .In a classification problem, f (x) = f (x ) means that x and x are not in the same class.The Lipschitz separation condition (1) becomes a margin condition specifying a minimum distance across classes: We can try to find a linear projection of x in some space V of lower dimension k, which separates f .It requires that f In most cases, the final dimension k can not be much smaller than d.
Linearization An alternative strategy is to linearize the variations of f with a first change of variable Φ(x) = {φ k (x)} k≤d of dimension d potentially much larger than the dimension d of x.We can then optimize a low-dimensional linear projection along directions where f is constant.We say that Φ separates f linearly if f (x) is well approximated by a one-dimensional projection: The regression vector w is optimized by minimizing a loss on the training data, which needs to be regularized if d > q, for example by an l p norm of w with a regularization constant λ: Sparse regressions are obtained with p ≤ 1, whereas p = 2 defines kernel regressions [16].
Classification problems are addressed similarly, by approximating the frontiers between classes.For example, a classification with Q classes can be reduced to Q − 1 "one versus all" binary classifications.Each binary classification is specified by an f (x) equal to 1 or −1 in each class.We approximate f (x) by f (x) = sign( Φ(x), w ), where w minimizes the training error (4).§3 Invariants, Symmetries and Diffeomorphisms We now study strategies to compute a change of variables Φ which linearizes f .Deep convolutional networks operate layer per layer and linearize f progressively, as depth increases.Classification and regression problems are addressed similarly by considering the level sets of f , defined by Ω t = {x : f (x) = t} if f is continuous.For classification, each level set is a particular class.Linear separability means that one can find w such that f (x) ≈ Φ(x), w .If x ∈ Ω t then Φ(x), w ≈ t, so all Ω t are mapped by Φ in different hyperplanes orthogonal to some w.The change of variable linearizes the level sets of f .Symmetries To linearize level sets, we need to find directions along which f (x) does not vary locally, and then linearize these directions in order to map them in a linear space.It is tempting to try to do this with some local data analysis along x.This is not possible because the training set includes few close neighbors in high dimension.We thus consider simultaneously all points x ∈ Ω and look for common directions along which f (x) does not vary.This is where groups of symmetries come in.Translations and diffeomorphisms will illustrate the difficulty to linearize high dimensional symmetries, and provide a first mathematical ground to analyze convolution networks architectures.
We look for invertible operators which preserve the value of f .The action of an operator g on x is written g.x.A global symmetry is an invertible and often non-linear operator g from Ω to Ω, such that f (g.x) = f (x) for all x ∈ Ω.If g 1 and g 2 are global symmetries then g 1 .g 2 is also a global symmetry, so products define groups of symmetries.Global symmetries are usually hard to find.We shall first concentrate on local symmetries.We suppose that there is a metric |g| G which measures the distance between g ∈ G and the identity.A function f is locally invariant to the action of G if We then say that G is a group of local symmetries of f .The constant C x is the local range of symmetries which preserve f .Since Ω is a continuous subset of R d , we consider groups of operators which transport vectors in Ω with a continuous parameter.They are called Lie groups if the group has a differential structure.
Translations and diffeomorphisms Let us interpolate the d samples of x and define x(u) for all u ∈ R n , with n = 1, 2, 3 respectively for time-series, images and volumetric data.The translation group G = R n is an example of Lie group.The action of g ∈ G = R n over x ∈ Ω is g.x(u) = x(u − g).The distance |g| G between g and the identity is the Euclidean norm of g ∈ R n .The function f is locally invariant to translations if sufficiently small translations of x do not change f (x).Deep convolutional networks compute convolutions, because they assume that translations are local symmetries of f .The dimension of a group G is the number of generators which define all group elements by products.For G = R n it is equal to n.
Translations are not powerful symmetries because they are defined by only n variables, and n = 2 for images.Many image classification problems are also locally invariant to small deformations, which provide much stronger constraints.It means that f is locally invariant to diffeomorphisms G = Diff(R n ), which transform x(u) with a differential warping of u ∈ R n .We do not know in advance what is the local range of diffeomorphism symmetries.For example, to classify images x of hand-written digits, certain deformations of x will preserve a digit class but modify the class of another digit.We shall linearize small diffeomorphims g.In a space where local symmetries are linearized, we can find global symmetries by optimizing linear projectors which preserve the values of f (x), and thus reduce dimensionality.
Local symmetries are linearized by finding a change of variable Φ(x) which locally linearizes the action of g ∈ G.We say that Φ is Lipschitz continuous if The norm x is just a normalization factor often set to 1.The Radon-Nikodim property proves that the map that transforms g into Φ(g.x) is almost everywhere differentiable in the sense of Gateaux.If |g| G is small then Φ(x) − Φ(g.x) is closely approximated by a bounded linear operator of g, which is the Gâteaux derivative.Locally, it thus nearly remains in a linear space.Lipschitz continuity over diffeomorphisms is defined relatively to a metric, which is now defined.A small diffeomorphism acting on x(u) can be written as a translation of u by a g(u): This diffeomorphism translates points by at most g ∞ = sup u∈R n |g(u)|.Let |∇g(u)| be the matrix norm of the Jacobian matrix of g at u. Small diffeomorphisms correspond to ∇g ).Their distance is thus multiplied by a scale factor, which is bounded above and below by 1 ± ∇g ∞ .The distance of this diffeomorphism to the identity is defined by: The factor 2 J is a local translation invariance scale.It gives the range of translations over which small diffeomorphisms are linearized.For J = ∞ the metric is globally invariant to translations.§4 Contractions and Scale Separation with Wavelets Deep convolutional networks can linearize the action of very complex non-linear transformations in high dimensions, such as inserting glasses in images of faces [28].A transformation of x ∈ Ω is a transport of x in Ω.To understand how to linearize any such transport, we shall begin with translations and diffeomorphisms.
Deep network architectures are covariant to translations, because all linear operators are implemented with convolutions.To compute invariants to translations and linearize diffeomorphisms, we need to separate scales and apply a non-linearity.This is implemented with a cascade of filters computing a wavelet transform, and a pointwise contractive non-linearity.Section 7 extends these tools to general group actions.
Averaging A linear operator can compute local invariants to the action of the translation group G, by averaging x along the orbit {g.x} g∈G , which are translations of x.This is done with a convolution by an averaging kernel φ One can verify [24] that this averaging is Lipschitz continuous to diffeomorphisms for all x ∈ L 2 (R n ), over a translation range 2 J .However, it eliminates the variations of x above the frequency 2 −J .If J = ∞ then Φ ∞ x = x(u) du, which eliminates nearly all information.
Wavelet transform A diffeomorphism acts as a local translation and scaling of the variable u.If we let aside translations for now, to linearize small diffeomorphism we must linearize this scaling action.This is done by separating the variations of x at different scales with wavelets.We define K wavelets ψ k (u) for u ∈ R n .They are regular functions with a fast decay and a zero average ψ k (u) du = 0.These K wavelets are dilated by 2 j : ψ j,k (u) = 2 −jn ψ k (2 −j u).A wavelet transform computes the local average of x at a scale 2 J , and variations at scales 2 j ≥ 2 J with wavelet convolutions: The parameter u is sampled on a grid such that intermediate sample values can be recovered by linear interpolations.The wavelets ψ k are chosen so that W is a contractive and invertible operator, and in order to obtain a sparse representation.This means that x ψ j,k (u) is mostly zero besides few high amplitude coefficients corresponding to variations of x(u) which "match" ψ k at the scale 2 j .This sparsity plays an important role in non-linear contractions.
For audio signals, n = 1, sparse representations are usually obtained with at least K = 12 intermediate frequencies within each octave 2 j , which are similar to half-tone musical notes.This is done by choosing a wavelet ψ(u) having a frequency bandwidth of less than 1/12 octave and For images, n = 2, we must discriminate image variations along different spatial orientation.It is obtained by separating angles πk/K, with an oriented wavelet which is rotated ψ k (u) = ψ(r −1 k u).Intermediate rotated wavelets are approximated by linear interpolations of these K wavelets.Figure 1 shows the wavelet transform of an image, with J = 4 scales and K = 4 angles, where x ψ j,k (u) is subsampled at intervals 2 j .It has few large amplitude coefficients shown in white.
Filter bank Wavelet transforms can be computed with a fast multiscale cascade of filters, which is at the core of deep network architectures.At each scale 2 j , we define a low-pass filter w j,0 which increases the averaging scale from 2 j−1 to 2 j , and band-pass filters w j,k which compute each wavelet: Let us write x j (u, 0) = x φ j (u) and x j (u, k) = x ψ j,k (u) for k = 0.It results from ( 11) that for 0 < j ≤ J and all 1 ≤ k ≤ K: These convolutions may be subsampled by 2 along u, in which case x j (u, k) is sampled at intervals 2 j along u.
Phase removal Wavelet coefficients x j (u, k) = x ψ j,k (u) oscillate at a scale 2 j .Translations of x smaller than 2 j modifies the complex phase of x j (u, k) if the wavelet is complex or its sign if it is real.Because of these oscillations, averaging x j with φ J outputs a zero signal.It is necessary to apply a non-linearity which removes oscillations.A modulus ρ(α) = |α| computes such a positive envelop.Averaging ρ(x ψ j,k (u)) by φ J outputs non-zero coefficients which are locally invariant at a scale 2 J : Replacing the modulus by a rectifier ρ(α) = max(0, α) gives nearly the same result, up to a factor 2. One can prove [24] that this representation is Lipschitz continuous to actions of diffeomorphisms over x ∈ L 2 (R n ), and thus satisfies (6) for the metric (8).Indeed, the wavelet coefficients of x deformed by g can be written as the wavelet coefficients of x with deformed wavelets.Small deformations produce small modifications of wavelets in L 2 (R n ), because they are localized and regular.The resulting modifications of wavelet coefficients is of the order of the diffeomorphism metric |g| Diff .
Figure 2: A convolution network iteratively computes each layer x j by transforming the previous layer x j−1 , with a linear operator W j and a pointwise non-linearity ρ.
Contractions A modulus and a rectifier are contractive non-linear pointwise operators: However, if α = 0 or α = 0 then this inequality is an equality.Replacing α and α by x ψ j,k (u) and x ψ j,k (u) shows that distances are much less reduced if x ψ j,k (u) is sparse.Such contractions do not reduce as much the distance between sparse signals and other signals.This is illustrated by reconstruction examples in Section 6.

Scale separation limitations
The local multiscale invariants in ( 13) have dominated pattern classification applications for music, speech and images, until 2010.It is called Mel-spectrum for audio [25] and SIFT type feature vectors [23] in images.Their limitations comes from the loss of information produced by the averaging by φ J in (13).To reduce this loss, they are computed at short time scales 2 J ≤ 50ms in audio signals, or over small image patches 2 2J = 16 2 pixels.As a consequence, they do not capture large scale structures, which are important for classification and regression problems.To build a rich set of local invariants at a large scale 2 J , it is not sufficient to separate scales with wavelets, we must also capture scale interactions.
A similar issue appears in physics to characterize the interactions of complex systems.Multiscale separations are used to reduce the parametrization of classical many body systems, for example with multipole methods [11].However, it does not apply to complex interactions, as in quantum systems.Interactions across scales, between small and larger structures, must be taken into account.Capturing these interactions with low-dimensional models is a major challenge.We shall see that deep neural networks and scattering transforms provide high order coefficients which partly characterize multiscale interactions.§5 Deep Convolutional Neural Network Architectures Deep convolutional networks are computational architectures introduced by Le Cun [20], providing remarkable regression and classification results in high dimension [21,19,17].We describe these architectures illustrated by Figure 2.They iterate over linear operators W j including convolutions, and predefined pointwise non-linearities.
A convolutional network takes in input a signal x(u), which is here an image.An internal network layer x j (u, k j ) at a depth j is indexed by the same translation variable u, usually subsampled, and a channel index k j .A layer x j is computed from x j−1 by applying a linear operator W j followed by a pointwise non-linearity ρ: x j = ρW j x j−1 .
The non-linearity ρ transforms each coefficient α of the array W j x j−1 , and satisfies the contraction condition (14).A usual choice is the rectifier ρ(α) = max(α, 0) for α ∈ R, but it can also be a sigmoid, or a modulus ρ(α) = |α| where α may be complex.
Since most classification and regression functions f (x) are invariant or covariant to translations, the architecture imposes that W j is covariant to translations.The output is translated if the input is translated.Since W j is linear, it can thus be written as a sum of convolutions: The variable u is usually subsampled.For a fixed j, all filters w j,kj (u, k) have the same support width along u, typically smaller than 10.
The operators ρ W j propagates the input signal x 0 = x until the last layer x J .This cascade of spatial convolutions defines translation covariant operators of progressively wider supports as the depth j increases.Each x j (u, k j ) is a non-linear function of x(v), for v in a square centered at u, whose width ∆ j does not depend upon k j .The width ∆ j is the spatial scale of a layer j.It is equal to 2 j ∆ if all filters w j,kj have a width ∆ and the convolutions (15) are subsampled by 2.
Neural networks include many side tricks.They sometimes normalize the amplitude of x j (v, k), by dividing it by the norm of all coefficients x j (v, k) for v in a neighborhood of u.This eliminates multiplicative amplitude variabilities.Instead of subsampling (15) on a regular grid, a max pooling may select the largest coefficients over each sampling cell.Coefficients may also be modified by subtracting a constant adapted to each coefficient.When applying a rectifier ρ, this constant acts as a soft threshold, which increases sparsity.It is usually observed that inside network coefficients x j (u, k j ) have a sparse activation.
The deep network output x J = Φ J (x) is provided to a classifier, usually composed of fully connected neural network layers [21].Supervised deep learning algorithms optimize the filter values w j,kj (u, k) in order to minimize the average classification or regression error on the training samples {x i , f (x i )} i≤q .There can be more than 10 8 variables in a network [21].The filter update is done with a back-propagation algorithm, which may be computed with a stochastic gradient descent, with regularization procedures such as dropout.This high-dimensional optimization is non-convex, but despite the presence of many local minima, the regularized stochastic gradient descent converges to a local minimum providing good accuracy on test examples [12].The rectifier non-linearity ρ is usually preferred because the resulting optimization has a better convergence.It however requires a large number of training examples.Several hundreds of examples per class are usually needed to reach a good performance.
Instabilities have been observed in some network architectures [31], where additions of small perturbations on x can produce large variations of x J .It happens when the norms of the matrices W j are larger than 1, and hence amplified when cascaded.However, deep network also have a strong form of stability illustrated by transfer learning [21].A deep network layer x J optimized on particular training databasis, can approximate different classification functions, if the final classification layers are trained on a new databasis.This means that it has learned stable structures, which can be transferred across similar learning problems.§6 Scattering on the Translation Group A deep network alternates linear operators W j and contractive non-linearities ρ.To analyze the properties of this cascade, we begin with a simpler architecture, where W j does not combine multiple convolutions across channels in each layer.We show that such network coefficients are obtained through convolutions with a reduced number of equivalent wavelet filters.It defines a scattering transform [24] whose contraction and linearization properties are reviewed.Variance reduction and loss of information are studied with reconstructions of stationary processes.
No channel combination Suppose that x j (u, k j ) is computed by convolving a single channel x j−1 (u, k j−1 ) along u: x j (u, k j ) = ρ x j−1 (•, k j−1 ) w j,h (u) with k j = (k j−1 , h) .
It corresponds to a deep network filtering (15), where filters do not combine several channels.Iterating on j defines a convolution tree, as opposed to a full network.It results from ( 16) that If ρ is a rectifier ρ(α) = max(α, 0) or a modulus ρ(α) = |α| then ρ(α) = α if α ≥ 0. We can thus remove this non-linearity at the output of an averaging filter w j,h .Indeed this averaging filter is applied to positive coefficients and thus computes positive coefficients, which are not affected by ρ.On the contrary, if w j,h is a band-pass filter then the convolution with x j−1 (•, k j−1 ) has alternating signs or a complex phase which varies.The non-linearity ρ removes the sign or the phase, which has a strong contraction effect.
Equivalent wavelet filter Let m be the number of band-pass filters {w jn,hj n } 1≤n≤m in the convolution cascade (17).All other filters are thus low-pass filters.If we remove ρ after each low-pass filter, we get m equivalent band-pass filters: The cascade of J convolutions ( 17) is reduced to m convolutions with these equivalent filters x J (u, k J ) = ρ(ρ(...ρ(ρ(x ψ j1,k1 ) ψ j2,k2 )...
with 0 < j 1 < j 2 < ... < j m−1 < J.If the final filter w J,h J at the depth J is a low-pass filter then ψ J,k J = φ J is an equivalent low-pass filter.In this case, the last non-linearity ρ can also be removed, which gives The operator Φ J x = x J is a wavelet scattering transform, introduced in [24].Changing the network filters w j,h modifies the equivalent band-pass filters ψ j,k .As in the fast wavelet transform algorithm (12), if w j,h is a rotation of a dilated filter w j then ψ j,h is a dilation and rotation of a single mother wavelet ψ.

Scattering order
The order m = 1 coefficients x J (u, k J ) = ρ(x ψ j1,k1 ) φ J (u) are the wavelet coefficient computed in (13).The loss of information due to averaging is now compensated by higher order coefficient.For m = 2, ρ(ρ(x ψ j1,k1 ) ψ j2,k2 ) φ J are complementary invariants.They measure interactions between variations of x at a scale 2 j1 , within a distance 2 j2 , and along orientation or frequency bands defined by k 1 and k 2 .These are scale interaction coefficients, missing from first order coefficients.Because ρ is strongly contracting, order m coefficients have an amplitude which decrease quickly as m increases [24,32].For images and audio signals, the energy of scattering coefficients becomes negligible for m ≥ 3. Let us emphasize that the convolution network depth is J, whereas m is the number of effective non-linearity of an output coefficient.
Diffeomorphism continuity Section 4 explains that a wavelet transform defines representations which are Lipschitz continuous to actions of diffeomorphisms.Scattering coefficients up to the order m are computed by applying m wavelet transforms.One can prove [24] that it thus defines a representation which is Lipschitz continuous to the the action of diffeomorphisms.There exists C > 0 such that plus a Hessian term which is neglected.This result is proved in [24] for ρ(α) = |α|, but it remains valid for any contractive pointwise operator such as rectifiers ρ(α) = max(α, 0).It relies on commutation properties of wavelet transforms and diffeomorphisms.It shows that the action of small diffeomorphisms is linearized over scattering coefficients.
Classification Scattering vectors are restricted to coefficients of order m ≤ 2, because their amplitude is negligible beyond.A translation scattering Φ J x is well adapted to classification problems where the main source of intra-class variability are due to translations, to small deformations, or to ergodic stationary processes.For example, intra-class variabilities of hand-written digit images are essentially due to translations and deformations.On the MNIST digit data basis [6], applying a linear classifier to scattering coefficients Φ J x gives state of the art classification errors.Music or speech classification over short time intervals of 100ms can be modeled by ergodic stationary processes.Good music and speech classification results are then obtained with a scattering transform [2].Image texture classification are also problems where intra class variability can be modeled by ergodic stationary processes.Scattering transforms give state of the art results over a wide range of image texture databases [6,29], compared to other descriptors including power spectrum moments.Softwares can be retrieved at www.di.ens.fr/data/software.

Stationary processes
To analyze the information loss, we now study the reconstruction of x from its scattering coefficients, in a stochastic framework where x is a stationary process.This will raise variance and separation issues, where sparsity plays a role.It also demonstrates the importance of second order scale interaction terms, to capture non-Gaussian geometric properties of ergodic stationary processes.Let us consider scattering coefficients of order m with φ J (u)du = 1.If x is a stationary process then ρ(...ρ(x ψ j1,k1 )... ψ jm,km ) remains stationary because convolutions and pointwise operators preserve stationarity.The spatial averaging by φ J provides a non-biased estimator of the expected value of Φ J x(u, k), which is a scattering moment: If x is a slow mixing process, which is a weak ergodicity assumption, then the estimation variance σ 2 J = Φ J x−E(Φ J x) 2 converges to zero [8] when J goes to ∞.Indeed, Φ J is computed by iterating on contractive operators, which average an ergodic stationary process x over progressively larger scales.One can prove that scattering moments characterize complex multiscale properties of fractals and multifractal processes, such as Brownian motions, Levi processes or Mandelbrot cascades [7].
Inverse scattering and sparsity Scattering transforms are generally not invertible but given Φ J (x) one can can compute vectors x such that Φ J (x) − Φ J (x) ≤ σ J .We initialize x0 as a Gaussian white noise realization, and iteratively update xn by reducing Φ J (x) − Φ J (x n ) with a gradient descent algorithm, until it reaches σ J [8].Since Φ J (x) is not convex, there is no guaranteed convergence, but numerical reconstructions converge up to a sufficient precision.The recovered x is a stationary process having nearly the same scattering moments as x, whose properties are similar to a maximum entropy process for fixed scattering moments [8].
Figure 3 shows several examples of images x with N 2 pixels.The first three images are realizations of ergodic stationary textures.The second row gives realizations of stationary Gaussian processes having the same N 2 second order covariance moments as the top textures.The third column shows the vorticity field of a two-dimensional turbulent fluid.The Gaussian realization is thus a Kolmogorov type model, which does not restore the filament geometry.The third row gives reconstructions from scattering coefficients, limited to order m ≤ 2. The scattering vector is computed at the maximum scale 2 J = N , with wavelets having K = 8 different orientations.It is thus completely invariant to translations.The dimension of Φ J x is about (K log 2 N ) 2 /2 N 2 .Scattering moments restore better texture geometries than the Gaussian models obtained with N 2 covariance moments.This geometry is mostly captured by second order scattering coefficients, providing scale interaction terms.Indeed, first order scattering moments can only reconstruct images which are similar to realizations of Gaussian processes.First and second order scattering moments also provide good models of ergodic audio textures [8].
The fourth image has very sparse wavelet coefficients.In this case the image is nearly perfectly restored by its scattering coefficients, up to a random translation.The reconstruction is centered for comparison.Section 4 explains that if wavelet coefficients are sparse then a rectifier or an absolute value contractions ρ does not contract as much distances with other signals.
Inverting a scattering transform is a non-linear inverse problem, which requires to recover a lost phase information.Sparsity has an important role on such phase recovery problems [32].Translating randomly the last motorcycle image defines a non-ergodic stationary process, whose wavelet coefficients are not as sparse.As a result, the reconstruction from a random initialization is very different, and does not preserve patterns which are important for most classification tasks.This is not surprising since there is much less scattering coefficients than image pixels.If we reduce 2 J so that the number of scattering coefficients reaches the number of pixels then the reconstruction is of good quality, but there is little variance reduction.
Concentrating on the translation group is not so effective to reduce variance when the process is not translation ergodic.Applying wavelet filters can destroy important structures which are not sparse over wavelets.Next section addresses both issues.Impressive texture synthesis results have been obtained with deep convolutional networks trained on image data bases [14], but with much more output coefficients.Numerical reconstructions [13] also show that one can also recover complex patterns, such as birds, airplanes, cars, dogs, ships, if the network is trained to recognize the corresponding image classes.The network keeps some form of memory of important classification patterns.§7 Multiscale Hierarchical Convolutional Networks Scattering transforms on the translation group are restricted deep convolutional network architectures, which suffer from variance issues and loss of information.We shall explain why channel combinations provide the flexibility needed to avoid some of these limitations.We analyze a general class of convolutional network architectures by extending the tools previously introduced.Contractions and invariants to translations are replaced by contractions along groups of local symmetries adapted to f , which are defined by parallel transports in each network layer.The network is structured by factorizing groups of symmetries, as depth increases.It implies that all linear operators can be written as generalized convolutions across multiple channels.To preserve the classification margin, wavelets must also be replaced by adapted filter weights, which separate discriminative patterns in multiple network fibers.
Separation margin Network layers x j = ρW j x j−1 are computed with operators ρW j which contract and separate components of x j .We shall see that W j also needs to prepare x j for the next transformation W j+1 , so consecutive operators W j and W j+1 are strongly dependant.Each W j is a contractive linear operator, W j z ≤ z to reduce the space volume, and avoid instabilities when cascading such operators [31].A layer x j−1 must separate f so that we can write f (x) = f j−1 (x j−1 ) for some function f j−1 (z).To simplify explanations, we concentrate on classification, where separation is an > 0 margin condition: The next layer x j = ρW j x j−1 lives in a contracted space but it must also satisfy The operator W j computes a linear projection which preserves this margin condition, but the resulting dimension reduction is limited.We can further contract the space non-linearly with ρ.To preserve the margin, it must reduce distances along non-linear displacements which transform any x j−1 into an x j−1 which is in the same class.
Parallel transport Displacements which preserve classes are defined by local symmetries (5), which are transformations ḡ such that f j−1 (x j−1 ) = f j−1 (ḡ.x j−1 ).To define a local invariant to a group of transformations G, we must process the orbit {ḡ.x j−1 } g∈G .However, W j is applied to x j−1 not on the non-linear transformations ḡ.x j−1 of x j−1 .The key idea is that a deep network can proceed in two steps.Let us write x j (u, k j ) = x j (v) with v ∈ P j .First, ρW j computes an approximate mapping of such an orbit {ḡ.x j−1 } ḡ∈G into a parallel transport in P j , which moves coefficients of x j .Then W j+1 applied to x j is filtering the orbits of this parallel transport.A parallel transport is defined by operators g ∈ G j acting on v ∈ P j , and we write The operator W j is defined so that G j is a group of local symmetries: f j (g.x j ) = f j (x j ) for small |g| Gj .This is obtained if a transport of x j = W j x j−1 by g ∈ G j corresponds to the action of a local symmetry ḡ of f j−1 By definition f j (x j ) = f j−1 (x j−1 ) = f (x).Since f j−1 (ḡ.x j−1 ) = f j−1 (x j−1 ) it results from ( 25) that f j (g.x j ) = f j (x j ).
The index space P j is called a G j -principal fiber bundle in differential geometry [26], illustrated by Figure 4.The orbits of G j in P j are fibers, indexed by the equivalence classes B j = P j /G j .They are globally invariant to the action of G j , and play an important role to separate f .Each fiber is indexing a continuous Lie group, but it is sampled along G j at intervals such that values of x j can be interpolated in between.As in the translation case, these sampling intervals depend upon the local invariance of x j , which increases with j.
Hierarchical symmetries In a hierarchical convolution network, we further impose that local symmetry groups are growing with depth, and can be factorized: The hierarchy begins for j = 0 by the translation group G 0 = R n , which acts on x(u) through the spatial variable u ∈ R n .The condition (26) is not necessarily satisfied by general deep networks, besides j = 0 for translations.It is used by joint scattering transforms [29,3] and has been proposed for unsupervised convolution network learning [9].Proposition 1 proves that this hierarchical embedding implies that each W j is a convolution on G j−1 .
Proposition 1.The group embedding (26) implies that x j can be indexed by (g, h, b) ∈ G j−1 × H j × B j and there exists w j,h.b ∈ C Pj−1 such that where h.b transports b ∈ B j by h ∈ H j in P j .
This proposition proves that W j is a convolution along the fibers of G j−1 in P j−1 .Each w j,h.b is a transformation of an elementary filter w j,b by a group of local symmetries h ∈ H j so that f j (x j (g, h, b)) remains constant when x j is locally transported along h.We give below several examples of groups H j and filters w j,h.b .However, learning algorithms compute filters directly, with no prior knowledge on the group H j .The filters w j,h.b can be optimized so that variations of x j (g, h, b) along h captures a large variance of x j−1 within each class.Indeed, this variance is then reduced by the next ρW j+1 .The generators of H j can be interpreted as principal symmetry generators, by analogy with the principal directions of a PCA.

Generalized scattering
The scattering convolution along translations ( 16) is replaced in ( 27) by a convolution along G j−1 , which combines different layer channels.Results for translations can essentially be extended to the general case.If w j,h.b is an averaging filter then it computes positive coefficients, so the non-linearity ρ can be removed.If each filter w j,h.b has a support in a single fiber indexed by b, as in Figure 4, then B j−1 ⊂ B j .It defines a generalized scattering transform, which is a structured multiscale hierarchical convolutional network such that G j−1 H j = G j and B j−1 ⊂ B j .If j = 0 then G 0 = P 0 = R n so B 0 is reduced to 1 fiber.
As in the translation case, we need to linearize small deformations in Diff(G j−1 ), which include much more local symmetries than the low-dimensional group G j−1 .A small diffeomorphism g ∈ Diff(G j−1 ) is a non-parallel transport along the fibers of G j−1 in P j−1 , which is a perturbation of a parallel transport.It modifies distances between pairs of points in P j−1 by scaling factors.To linearize such diffeomorphisms, we must use localized filters whose supports have different scales.Scale parameters are typically different along the different generators of G j−1 = R n H 1 ... H j−1 .Filters can be constructed with wavelets dilated at different scales, along the generators of each group H k for 1 ≤ k ≤ j.Linear dimension reduction mostly results from this filtering.Variations at fine scales may be eliminated, so that x j (g, h, b) can be coarsely sampled along g.

Rigid movements
For small j, the local symmetry groups H j may be associated to linear or non-linear physical phenomena such as rotations, scaling, colored illuminations or pitch frequency shifts.Let SO(n) be the group of rotations.Rigid movements SE(n) = R n SO(n) is a non-commutative group, which often includes local symmetries.For images, n = 2, this group becomes a transport in P 1 with H 1 = SO(n) which rotates a wavelet filter w 1,h (u) = w 1 (r −1 h u).Such filters are often observed in the first layer of deep convolutional networks [13].They map the action of ḡ = (v, r k ) ∈ SE(n) on x to a parallel transport of (u, h) ∈ P 1 defined for g ∈ G 1 = R 2 ×SO(n) by g.(u, h) = (v +r k u, h+k).Small diffeomorphisms in Diff(G j ) correspond to deformations along translations and rotations, which are sources of local symmetries.A rototranslation scattering [29,27] linearizes them with wavelet filters along translations and rotations, with G j = SE(n) for all j > 1.This roto-translation scattering can efficiently regress physical functionals which are often invariant to rigid movements, and Lipschitz continuous to deformations.For example, quantum molecular energies f (x) are well estimated by sparse regressions over such scattering representations [18].
Audio pitch Pitch frequency shift is a more complex example of a non-linear symmetry for audio signals.Two different musical notes of a same instrument have a pitch shift.Their harmonic frequencies are multiplied by a factor 2 h , but it is not a dilation because the note duration is not changed.With narrow band-pass filters w 1,h (u) = w 1 (2 −h u), a pitch shift is approximatively mapped to a translation along h ∈ H 1 = R of ρ(x w 1,h (u)), with no modification along the time u.The action of g A pitch shift also comes with deformations along time and log-frequencies, which define a much larger class of symmetries in Diff(G 1 ).Two-dimensional wavelets along (u, h) can linearize these small time and log-frequency deformations.These define a joint time-frequency scattering applied to speech and music classifications [3].Such transformations were first proposed as neurophysiological models of audition [25].

Manifolds of patterns
The group H j is associated to complex transformations when j increases.It needs to capture large transformations between different patterns in a same class, for example chairs of different styles.Let us consider training samples {x i } i of a same class.The iterated network contractions transform them into vectors {x i j−1 } i which are much closer.Their distances define weighted graphs which sample underlying continuous manifolds in the space.Such manifolds clearly appear in [5], for high-level patterns such as chairs or cars, together with poses and colors.As opposed to manifold learning, deep network filters result from a global optimization which can be computed in high dimension.The principal symmetry generators of H j is associated to common transformations over all manifolds of examples x i j−1 , which preserve the class while capturing large intra-class variance.They are approximatively mapped to a parallel transport in x j by the filters w j,h.b .The diffeomorphisms in Diff(G j ) are non-parallel transports corresponding to high-dimensional displacements on the manifolds of x j−1 .Linearizing Diff(G j ) is equivalent to partially flatten simultaneously all these manifolds, which may explain why manifolds are progressively more regular as the network depth increases [5], but it involves open mathematical questions.
Sparse support vectors We have up to now been concentrated on the reduction of the data variability through contractions.We now explain why the classification margin can be preserved thanks to the existence of multiple fibers B j in P j , by adapting filters instead of using standard wavelets.The fibers indexed by b ∈ B j are separation instruments, which increase dimensionality to avoid reducing the classification margin.They prevent from collapsing vectors in different classes, which have a distance x j−1 − x j−1 close to the minimum margin .These vectors are close to classification frontiers.They are called multiscale support vectors, by analogy with support vector machines.To avoid further contracting their distance, they can be separated along different fibers indexed by b.The separation is achieved by filters w j,h.b , which transform x j−1 and x j−1 into x j (g, h, b) and x j (g, h, b) having sparse supports on different fibers b.The next contraction ρW j+1 reduces distances along fibers indexed by (g, h) ∈ G j , but not across b ∈ B j , which preserves distances.The contraction increases with j so the number of support vectors close to frontiers also increases, which implies that more fibers are needed to separate them.When j increases, the size of x j is a balance between the dimension reduction along fibers, by subsampling g ∈ G j , and an increasing number of fibers B j which encode progressively more support vectors.Coefficients in these fibers become more specialized and invariants, as the grandmother neurons observed in deep layers of convolutional networks [1].They have a strong response to particular patterns and are invariant to a large class of transformations.In this model, the choice of filters w j,h.b are adapted to produce sparse representations of multiscale support vectors.They provide a sparse distributed code, defining invariant pattern memorisation.This memorisation is numerically observed in deep network reconstructions [13], which can restore complex patterns within each class.Let us emphasize that groups and fibers are mathematical ghost behind filters, which are never computed.The learning optimization is directly performed on filters, which carry the trade-off between contractions to reduce the data variability and separation to preserve classification margin.§8 Conclusion This paper provides a mathematical framework to analyze contraction and separation properties of deep convolutional networks.In this model, network filters are guiding non-linear contractions, to reduce the data variability in directions of local symmetries.The classification margin can be controlled by sparse separations along network fibers.Network fibers combine invariances along groups of symmetries and distributed pattern representations, which could be sufficiently stable to explain transfer learning of deep networks [21].However, this is only a framework.We need complexity measures, approximation theorems in spaces of high-dimensional functions, and guaranteed convergence of filter optimization, to fully understand the mathematics of these convolution networks.
Besides learning, there are striking similarities between these multiscale mathematical tools and the treatment of symmetries in particle and statistical physics [15].One can expect a rich cross fertilization between high-dimensional learning and physics, through the development of a common mathematical language.

Figure 1 :
Figure 1: Wavelet transform of an image x(u), computed with a cascade of convolutions with filters over J = 4 scales and K = 4 orientations.The low-pass and K = 4 band-pass filters are shown on the first arrows.

Figure 3 :
Figure 3: First row: original images.Second row: realization of a Gaussian process with same second covariance moments.Third row: reconstructions from first and second order scattering coefficients.

Figure 4 :
Figure 4: A multiscale hierarchical networks computes convolutions along the fibers of a parallel transport.It is defined by a group G j of symmetries acting on the index set P j of a layer x j .Filter weights are transported along fibers.