## Abstract

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 [1]. They give state-of-the-art performances for image classification with thousands of complex classes [2], speech recognition [3], biomedical applications [4], natural language understanding [5] and in many other domains. They are also studied as neurophysiological models of vision [6].

Multilayer neural networks are computational learning architectures that propagate the input data across a sequence of linear operators and simple nonlinearities. The properties of shallow networks, with one hidden layer, are well understood as decompositions in families of ridge functions [7]. However, these approaches do not extend to networks with more layers. Deep convolutional neural networks, introduced by Le Cun [8], are implemented with linear convolutions followed by nonlinearities, over typically more than five 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 [1,6], but relations with networks weights and nonlinearities are complex. This paper aims at clarifying important principles that govern the properties of such networks, although 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 one 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 problems.

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 §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 §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 §7.

To linearize diffeomorphisms and preserve separability, §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 §5. They iterate on linear operators which filter and linearly combine different channels in each network layer, followed by contractive nonlinearities.

To understand how nonlinear contractions interact with linear operators, §6 begins with simpler networks which do not recombine channels in each layer. It defines a nonlinear scattering transform, introduced in [9], 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 fibres which are progressively more invariant and specialized. This can give rise to invariant grandmother type neurons observed in deep networks [10]. The paper studies architectures as opposed to computational learning of network weights, which is an outstanding optimization issue [1].

*Notations*. ∥*z*∥ is a Euclidean norm if *z* is a vector in a Euclidean space. If *z* is a function in **L**^{2}, then . If *z*={*z*_{k}}_{k} is a sequence of vectors or functions, then .

## 2. Linearization, projection and separability

Supervised learning computes an approximation 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 , not a low-dimensional manifold. In a regression problem, *f*(*x*) takes its values in , 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 *f*(*x*)=*f*_{0}(*Φ*(*x*)). It is equivalent to impose that if *f*(*x*)≠*f*(*x*′) then *Φ*(*x*)≠*Φ*(*x*′). We then say that *Φ* *separates* *f*. For regression problems, to guarantee that *f*_{0} is regular, we further impose that the separation is Lipschitz:
2.1It implies that *f*_{0} is Lipschitz continuous: |*f*_{0}(*z*)−*f*_{0}(*z*′)|≤*ϵ*^{−1}|*z*−*z*′|, for (*z*,*z*′)∈*Φ*(*Ω*)^{2}. In a classification problem, *f*(*x*)≠*f*(*x*′) means that *x* and *x*′ are not in the same class. The Lipschitz separation condition (2.1) becomes a margin condition specifying a minimum distance across classes:
2.2We can try to find a linear projection of *x* in some space **V** of lower dimension *k*, which separates *f*. It requires that *f*(*x*)=*f*(*x*+*z*) for all *z*∈**V**^{⊥}, where **V**^{⊥} is the orthogonal complement of **V** in , of dimension *d*−*k*. In most cases, the final dimension *k* cannot 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: 2.3The 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

*λ*: 2.4Sparse regressions are obtained with

*p*≤1, whereas

*p*=2 defines kernel regressions [11].

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 , where *w* minimizes the training error (2.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*〉. It implies that *Φ*(*Ω*_{t}) is in a hyperplane orthogonal to some *w*.

*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 neighbours 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 analyse convolution networks architectures.

We look for invertible operators that 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 nonlinear 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
3.1We 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*. Because *Ω* is a continuous subset of , 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 , with *n*=1,2,3, respectively, for time-series, images and volumetric data. The translation group is an example of a Lie group. The action of over *x*∈*Ω* is *g*.*x*(*u*)=*x*(*u*−*g*). The distance |*g*|_{G} between *g* and the identity is the Euclidean norm of . 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 that define all group elements by products. For 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 , which transform *x*(*u*) with a differential warping of . We do not know in advance what is the local range of diffeomorphism symmetries. For example, to classify images *x* of handwritten 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 that 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
3.2The 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 Gâteaux. 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 relative 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*):
3.3This diffeomorphism translates points by at most . Let |∇*g*(*u*)| be the matrix norm of the Jacobian matrix of *g* at *u*. Small diffeomorphisms correspond to . Applying a diffeomorphism *g* transforms two points (*u*_{1},*u*_{2}) into (*u*_{1}−*g*(*u*_{1}),*u*_{2}−*g*(*u*_{2})). Their distance is thus multiplied by a scale factor, which is bounded above and below by . The distance of this diffeomorphism to the identity is defined by
3.4The factor 2^{J} is a local translation invariance scale. It gives the range of translations over which small diffeomorphisms are linearized. For , the metric is globally invariant to translations.

## 4. Contractions and scale separation with wavelets

Deep convolutional networks can linearize the action of very complex nonlinear transformations in high dimensions, such as inserting glasses in images of faces [12]. 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 nonlinearity. This is implemented with a cascade of filters computing a wavelet transform, and a pointwise contractive nonlinearity. Section 7 extends these tools to general group actions.

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 *ϕ*_{J}(*u*)=2^{−nJ}*ϕ*(2^{−J}*u*) of size 2^{J}, with :
4.1One can verify [9] that this averaging is Lipschitz continuous to diffeomorphisms for all , over a translation range 2^{J}. However, it eliminates the variations of *x* above the frequency 2^{−J}. If then , 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 a small diffeomorphism, then 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 . They are regular functions with a fast decay and a zero average . 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:
4.2The 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 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 nonlinear 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 *ψ*_{k}(*u*)=2^{k/K}*ψ*(2^{−k/K}*u*) for 1≤*k*≤*K*. 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 . 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:
4.3Let us write *x*_{j}(*u*,0)=*x*★*ϕ*_{j}(*u*) and *x*_{j}(*u*,*k*)=*x*★*ψ*_{j,k}(*u*) for *k*≠0. It results from (4.3) that for 0<*j*≤*J* and all 1≤*k*≤*K*:
4.4These convolutions may be subsampled by 2 along *u*; in that case *x*_{j}(*u*,*k*) is sampled at intervals 2^{j} along *u*.

*Phase removal contraction.* Wavelet coefficients *x*_{j}(*u*,*k*)=*x*★*ψ*_{j,k}(*u*) oscillate at a scale 2^{j}. Translations of *x* smaller than 2^{j} modify 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 nonlinearity which removes oscillations. A modulus *ρ*(*α*)=|*α*| computes such a positive envelope. Averaging *ρ*(*x*★*ψ*_{j,k}(*u*)) by *ϕ*_{J} outputs non-zero coefficients which are locally invariant at a scale 2^{J}:
4.5Replacing the modulus by a rectifier gives nearly the same result, up to a factor 2. One can prove [9] that this representation is Lipschitz continuous to actions of diffeomorphisms over , and thus satisfies (3.2) for the metric (3.4). 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 , because they are localized and regular. The resulting modifications of wavelet coefficients are of the order of the diffeomorphism metric |*g*|_{Diff}.

A modulus and a rectifier are contractive nonlinear pointwise operators:
4.6However, 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 §6.

*Scale separation limitations.* The local multiscale invariants in (4.5) have dominated pattern classification applications for music, speech and images, until 2010. It is called *Mel-spectrum* for audio [13] and SIFT-type feature vectors [14] in images. Their limitations comes from the loss of information produced by the averaging by *ϕ*_{J} in (4.5). To reduce this loss, they are computed at short time scales 2^{J}≤50 ms 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 [15]. 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 [8], providing remarkable regression and classification results in high dimension [1–3]. We describe these architectures illustrated by figure 2. They iterate over linear operators *W*_{j} including convolutions, and predefined pointwise nonlinearities.

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 nonlinearity *ρ*:
The nonlinearity *ρ* transforms each coefficient *α* of the array *W*_{j}*x*_{j−1}, and satisfies the contraction condition (4.6). A usual choice is the rectifier for , but it can also be a sigmoid, or a modulus *ρ*(*α*)=|*α*| where *α* may be complex.

Because 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. Because *W*_{j} is linear, it can thus be written as a sum of convolutions:
5.1The 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} propagate 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 nonlinear function of *x*(*v*), for *v* in a square centred 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 (5.1) 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 neighbourhood of *u*. This eliminates multiplicative amplitude variabilities. Instead of subsampling (5.1) 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 [1]. 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 [1]. 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 [16]. The rectifier nonlinearity *ρ* 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 [17], 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 networks also have a strong form of stability illustrated by transfer learning [1]. A deep network layer *x*_{J}, optimized on particular training databases, can approximate different classification functions, if the final classification layers are trained on a new database. 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 nonlinearities *ρ*. To analyse 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 [9] 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*:
6.1It corresponds to a deep network filtering (5.1), where filters do not combine several channels. Iterating on *j* defines a convolution tree, as opposed to a full network. It results from (6.1) that
6.2If *ρ* is a rectifier or a modulus *ρ*(*α*)=|*α*| then *ρ*(*α*)=*α* if *α*≥0. We can thus remove this nonlinearity 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 nonlinearity *ρ* removes the sign or the phase, which has a strong contraction effect.

*Equivalent wavelet filter.* Suppose that there are *m* band-pass filters {*w*_{jn,hjn}}_{1≤n≤m} in the convolution cascade (6.2), and that all others are low-pass filters. If we remove *ρ* after each low-pass filter, then we get *m* equivalent band-pass filters:
6.3The cascade of *J* convolutions (6.2) is reduced to *m* convolutions with these equivalent filters
6.4with 0<*j*_{1}<*j*_{2}<⋯<*j*_{m−1}<*J*. If the final filter *w*_{J,hJ} at the depth *J* is a low-pass filter, then *ψ*_{J,kJ}=*ϕ*_{J} is an equivalent low-pass filter. In this case, the last nonlinearity *ρ* can also be removed, which gives
6.5The operator *Φ*_{J}*x*=*x*_{J} is a wavelet scattering transform, introduced in [9]. Changing the network filters *w*_{j,h} modifies the equivalent band-pass filters *ψ*_{j,k}. As in the fast wavelet transform algorithm (4.4), 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 (4.5). The loss of information owing 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 decreases quickly as *m* increases [9,18]. 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 nonlinearity of an output coefficient.

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 [9] that it thus defines a representation which is Lipschitz continuous to the action of diffeomorphisms. There exists *C*>0 such that
plus a Hessian term which is neglected. This result is proved in [9] for *ρ*(*α*)=|*α*|, but it remains valid for any contractive pointwise operator such as rectifiers . 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 sources of intraclass variability are owing to translations, to small deformations or to ergodic stationary processes. For example, intraclass variabilities of handwritten digit images are essentially owing to translations and deformations. On the MNIST digit data basis [19], 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 100 ms can be modelled by ergodic stationary processes. Good music and speech classification results are then obtained with a scattering transform [20]. Image texture classification is also problems where intraclass variability can be modelled by ergodic stationary processes. Scattering transforms give state-of-the-art results over a wide range of image texture databases [19,21], compared with other descriptors including power spectrum moments.

*Stationary processes.* To analyse 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*
6.6with . 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:
6.7If *x* is a slow mixing process, which is a weak ergodicity assumption, then the estimation variance converges to zero [22] 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 [23].

*Inverse scattering and sparsity.* Scattering transforms are generally not invertible but given *Φ*_{J}(*x*) one can compute vectors such that . We initialize as a Gaussian white noise realization, and iteratively update by reducing with a gradient descent algorithm, until it reaches *σ*_{J} [22]. Because *Φ*_{J}(*x*) is not convex, there is no guaranteed convergence, but numerical reconstructions converge up to a sufficient precision. The recovered is a stationary process having nearly the same scattering moments as *x*, the properties of which are similar to a maximum entropy process for fixed scattering moments [22].

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 . 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 [22].

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 centred for comparison. Section 4 explains that if wavelet coefficients are sparse then a rectifier or an absolute value contractions *ρ* does not contract as much as distances with other signals. Indeed, |*ρ*(*α*)−*ρ*(*α*′)|=|*α*−*α*′| if *α*=0 or *α*′=0. Inverting a scattering transform is a nonlinear inverse problem, which requires one to recover a lost phase information. Sparsity has an important role in such phase recovery problems [18]. 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, because there are much fewer 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. Section 7 addresses both issues. Impressive texture synthesis results have been obtained with deep convolutional networks trained on image databases [24], but with much more output coefficients. Numerical reconstructions [25] also show that one can also recover complex patterns, such as birds, aeroplanes, 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 analyse 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 fibres.

*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 dependent. Each *W*_{j} is a contractive linear operator, ∥*W*_{j}*z*∥≤∥*z*∥ to reduce the space volume, and avoid instabilities when cascading such operators [17]. 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:
7.1The next layer *x*_{j}=*ρW*_{j}*x*_{j−1} lives in a contracted space but it must also satisfy
7.2The 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 nonlinearly with *ρ*. To preserve the margin, it must reduce distances along nonlinear displacements which transform any *x*_{j−1} into an *x*′_{j−1} which is in the same class. Such displacements are defined by local symmetries (3.1), which are transformations such that .

*Parallel transport.* To define a local invariant to a group of transformations *G*, we must process the orbit . However, *W*_{j} is applied to *x*_{j−1} not on the nonlinear transformations 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 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 transports by *g*∈*G*_{j} correspond approximatively to local symmetries of *f*_{j−1}. Indeed, if then .

The index space *P*_{j} is called a *G*_{j}-principal fibre bundle in differential geometry [26], illustrated by figure 4. The orbits of *G*_{j} in *P*_{j} are fibres, 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 fibre 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:
7.3The hierarchy begins for *j*=0 by the translation group , which acts on *x*(*u*) through the spatial variable . The condition (7.3) is not necessarily satisfied by general deep networks, besides *j*=0 for translations. It is used by joint scattering transforms [21,27] and has been proposed for unsupervised convolution network learning [28]. Proposition 7.1 proves that this hierarchical embedding implies that each *W*_{j} is a convolution on *G*_{j−1}.

### Proposition 7.1

*The group embedding (7.3) implies that x*_{j} *can be indexed by (g,h,b)∈G*_{j−1}*×H*_{j}*×B*_{j} *and there exists* *such that
*7.4*where h.b transports b∈B*_{j} *by h∈H*_{j} *in P*_{j}.

### Proof.

We write *x*_{j}=*ρW*_{j}*x*_{j−1} as inner products with row vectors:
7.5If then . One can write with and *b*∈*B*_{j}=*P*_{j}/*G*_{j}. If then can be decomposed into , where *g*.*x*_{j}=*ρ*(〈*g*.*x*_{j−1},*w*_{j,b}〉). But *g*.*x*_{j−1}(*v*′)=*x*_{j−1}(*g*.*v*′) so with a change of variable, we get *w*_{j,g.b}(*v*′)=*w*_{j,b}(*g*^{−1}.*v*′). Hence, . Inserting this filter expression in (7.5) proves (7.4). ▪

This proposition proves that *W*_{j} is a convolution along the fibres 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* capture 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 principal component analysis.

*Generalized scattering.* The scattering convolution along translations (6.1) is replaced in (7.4) 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 nonlinearity *ρ* can be removed. If each filter *w*_{j,h.b} has a support in a single fibre 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 and *B*_{j−1}⊂*B*_{j}. If *j*=0, then so *B*_{0} is reduced to 1 fibre.

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 fibres 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 . 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 with linear or nonlinear physical phenomena such as rotations, scaling, coloured illuminations or pitch frequency shifts. Let *SO*(*n*) be the group of rotations. Rigid movements 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 . Such filters are often observed in the first layer of deep convolutional networks [25]. They map the action of on *x* to a parallel transport of (*u*,*h*)∈*P*_{1} defined for 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 roto-translation scattering [21,29] 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 [30].

*Audio pitch.* Pitch frequency shift is a more complex example of a nonlinear 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 of *ρ*(*x*★*w*_{1,h}(*u*)), with no modification along the time *u*. The action of over (*u*,*h*)∈*P*_{1} is thus a two-dimensional translation *g*.(*u*,*h*)=(*u*+*v*,*h*+*k*). 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 [27]. Such transformations were first proposed as neurophysiological models of audition [13].

*Manifolds of patterns.* The group *H*_{j} is associated with 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 which are much closer. Their distances define weighted graphs which sample underlying continuous manifolds in the space. Such manifolds clearly appear in [31], for high-level patterns such as chairs or cars, together with poses and colours. 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} are associated with common transformations over all manifolds of examples , which preserve the class while capturing large intraclass 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 flattening simultaneously all these manifolds, which may explain why manifolds are progressively more regular as the network depth increases [31], but it involves open mathematical questions.

*Sparse support vectors.* We have up to now 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 fibres *B*_{j} in *P*_{j}, by adapting filters instead of using standard wavelets. The fibres 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 fibres 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 fibres *b*. The next contraction *ρW*_{j+1} reduces distances along fibres 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 fibres are needed to separate them.

When *j* increases, the size of *x*_{j} is a balance between the dimension reduction along fibres, by subsampling *g*∈*G*_{j}, and an increasing number of fibres *B*_{j} which encode progressively more support vectors. Coefficients in these fibres become more specialized and invariants, as the grandmother neurons observed in deep layers of convolutional networks [10]. 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 memorization. This memorization is numerically observed in deep network reconstructions [25], which can restore complex patterns within each class. Let us emphasize that groups and fibres are mathematical ghosts 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 analyse contraction and separation properties of deep convolutional networks. In this model, network filters are guiding nonlinear contractions, to reduce the data variability in directions of local symmetries. The classification margin can be controlled by sparse separations along network fibres. Network fibres combine invariances along groups of symmetries and distributed pattern representations, which could be sufficiently stable to explain transfer learning of deep networks [1]. 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 [32]. One can expect a rich cross fertilization between high-dimensional learning and physics, through the development of a common mathematical language.

## Data accessibility

Softwares reproducing experiments can be retrieved at www.di.ens.fr/data/software.

## Competing interests

The author declares that there are no competing interests.

## Funding

This work was supported by the ERC grant InvariantClass 320959.

## Acknowledgements

I thank Carmine Emanuele Cella, Ivan Dokmaninc, Sira Ferradans, Edouard Oyallon and Irène Waldspurger for their helpful comments and suggestions.

## Footnotes

One contribution of 13 to a theme issue ‘Adaptive data analysis: theory and applications’.

- Accepted December 16, 2015.

- © 2016 The Authors.