Invariant Scattering Convolution Networks

Introduction

A major difficulty of image classification comes from the considerable variability within image classes and the inability of Euclidean distances to measure image similarities. Part of this variability is due to rigid translations, rotations or scaling. This variability is often uninformative for classification and should thus be eliminated. In the framework of kernel classifiers , metrics are defined as a Euclidean distance applied on a representation $\Phi(x)$ of signals $x$. The operator $\Phi$ must therefore be invariant to these rigid transformations.

Non-rigid deformations also induce important variability within object classes . For instance, in handwritten digit recognition, one must take into account digit deformations due to different writing styles. However, a full deformation invariance would reduce discrimination since a digit can be deformed into a different digit, for example a one into a seven. The representation must therefore not be deformation invariant but continuous to deformations, to handle small deformations with a kernel classifier. A small deformation of an image $x$ into $x^{\prime}$ should correspond to a small Euclidean distance $\|\Phi(x)-\Phi(x^{\prime})\|$ in the representation space, as further explained in Section 2.

Translation invariant representations can be constructed with registration algorithms or with the Fourier transform modulus. However, Section 2.1 explains why these invariants are not stable to deformations and hence not adapted to image classification. Trying to avoid Fourier transform instabilities suggests replacing sinusoidal waves by localized waveforms such as wavelets. However, wavelet transforms are not invariant to translations. Building invariant representations from wavelet coefficients requires introducing non-linear operators, which leads to a convolution network architecture.

Deep convolution networks have the ability to build large-scale invariants which are stable to deformations . They have been applied to a wide range of image classification tasks. Despite the remarkable successes of this neural network architecture, the properties and optimal configurations of these networks are not well understood because of cascaded non-linearities. Why use multiple layers ? How many layers ? How to optimize filters and pooling non-linearities ? How many internal and output neurons ? These questions are mostly answered through numerical experimentations that require significant expertise.

Deformation stability is obtained with localized wavelet filters which separate the image variations at multiple scales and orientations . Computing a non-zero translation invariant representation from wavelet coefficients requires introducing a non-linearity, which is chosen to be a modulus to optimize stability . Wavelet scattering networks, introduced in , build translation invariant representations with average poolings of wavelet modulus coefficients. The output of the first network layer is similar to SIFT or Daisy type descriptors. However, this limited set of locally invariant coefficients is not sufficiently informative to discriminate complex structures over large-size domains. The information lost by the averaging is recovered by computing a next layer of invariant coefficients, with the same wavelet convolutions and average modulus poolings. A wavelet scattering is thus a deep convolution network which cascades wavelet transforms and modulus operators. The mathematical properties of scattering operators explain how these deep network coefficients relate to image sparsity and geometry. The network architecture is optimized in Section 3, to retain important information while avoiding useless computations.

A scattering representation of stationary processes is introduced for texture discrimination. As opposed to the Fourier power spectrum, it provides information on higher order moments and can thus discriminate non-Gaussian textures having the same power spectrum. Classification applications are studied in Section 4.1. Scattering classification properties are demonstrated with a Gaussian kernel SVM and a generative classifier, which selects affine space models computed with a PCA. State-of-the-art results are obtained for handwritten digit recognition on MNIST and USPS databes, and for texture discrimination. Software is available at www.cmap.polytechnique.fr/scattering.

Towards a Convolution Network

Section 2.1 formalizes the deformation stability condition as a Lipschitz continuity property, and explains why high Fourier frequencies are source of unstabilites. Section 2.2 introduces a wavelet-based scattering transform, which is translation invariant and stable to deformations, and section 2.3 describes its convolutional network architecture.

A canonical invariant $\Phi(x)=x(u-a(x))$ registers $x$ with an anchor point $a(x)$, which is translated when $x$ is translated: $a(L_{c}x)=a(x)+c$. It is therefore invariant: $\Phi(L_{c}x)=\Phi(x)$. For example, the anchor point may be a filtered maxima $a(x)=\arg\max_{u}|x\star h(u)|$, for some filter $h(u)$.

The Fourier transform modulus is another example of translation invariant representation. Let $\hat{x}(\omega)$ be the Fourier transform of $x(u)$. Since $\widehat{L_{c}x}(\omega)=e^{-ic.\omega}\,\hat{x}(\omega)$, it results that $|\widehat{L_{c}x}|=|\hat{x}|$ does not depend upon $c$.

To obtain appropriate similarity measurements between images which have undergone non-rigid transformations, the representation must also be stable to small deformations. A small deformation can be written $L_{\tau}x(u)=x(u-\tau(u))$ where $\tau(u)$ depends upon $u$ and thus deforms the image. The deformation gradient tensor $\nabla\tau(u)$ is a matrix whose norm $|\nabla\tau(u)|$ measures the deformation amplitude at $u$. A small deformation is an invertible transformation if $|\nabla\tau(u)|<1$ . Stability to deformations is expressed as a Lipschitz continuity condition relative to this deformation metric:

where $\|x\|^{2}=\int|x(u)|^{2}\,du$. This property implies global translation invariance, because if $\tau(u)=c$ then $\nabla\tau(u)=0$, but it is much stronger.

A Fourier modulus is translation invariant but unstable with respect to deformations at high frequencies. Indeed, $|\,|\hat{x}(\omega)|-|\widehat{L_{\tau}x}(\omega)|\,|$ can be arbitrarily large at a high frequency $\omega$, even for small deformations and in particular small dilations. As a result, $\Phi(x)=|\hat{x}|$ does not satisfy the deformation continuity condition (2) . A Fourier modulus also loses too much information. For example, a Dirac $\delta(u)$ and a linear chirp $e^{iu^{2}}$ are totally different signals having Fourier transforms whose moduli are equal and constant. Very different signals may not be discriminated from their Fourier modulus.

A registration invariant $\Phi(x)=x(u-a(x))$ carries more information than a Fourier modulus, and characterizes $x$ up to a global absolute position information . However, it has the same high-frequency instability as a Fourier transform. Indeed, for any choice of anchor point $a(x)$, applying the Plancherel formula proves that

If $x^{\prime}=L_{\tau}x$, the Fourier transform instability at high frequencies implies that $\Phi(x)=x(u-a(x))$ is also unstable with respect to deformations.

2 Scattering Wavelets

A wavelet is a localized waveform and is thus stable to deformation, as opposed to the Fourier sinusoidal waves. A wavelet transform computes convolutions with wavelets. It is thus translation covariant, not invariant. A scattering transform computes non-linear invariants with modulus and averaging pooling functions.

A wavelet transform filters $x$ using a family of wavelets: $\{x\star\psi_{\lambda}(u)\}_{\lambda}$. It is computed with a filter bank of dilated and rotated wavelets having no orthogonality property. As further explained in Section 3.1, it is stable and invertible if the rotated and scaled wavelet filters cover the whole frequency plane. On discrete images, to avoid aliasing, we only capture frequencies in the circle $|\omega|\leq\pi$ inscribed in the image frequency square. However, most digital natural images and textures have negligible energy outside this frequency circle.

where $C_{2}$ is adjusted so that $\int\psi(u)\,du=0$. Figure 1 shows the Morlet wavelet with $\sigma=0.85$ and $\xi=3\pi/4$, used in all classification experiments.

A wavelet transform commutes with translations, and is therefore not translation invariant. To build a translation invariant representation, it is necessary to introduce a non-linearity. If $R$ is a linear or non-linear operator which commutes with translations, $R(L_{c}x)=L_{c}Rx$, then the integral $\int Rx(u)\,du$ is translation invariant. Applying this to $Rx=x\star\psi_{\lambda}$ gives a trivial invariant $\int x\star\psi_{\lambda}(u)\,du=0$ for all $x$ because $\int\psi_{\lambda}(u)\,du=0$. If $Rx=M(x\star\psi_{\lambda})$ but $M$ is linear and commutes with translations then the integral still vanishes, which imposes choosing a non-linear $M$. Taking advantage of the wavelet transform stability to deformations, to obtain integrals which are also stable to deformations we also impose that $M$ commutes with deformations

More translation invariant coefficients can be computed by further iterating on the wavelet transform and modulus operators. Let $U[\lambda]x=|x\star\psi_{\lambda}|$. Any sequence $p=(\lambda_{1},\lambda_{2},...,\lambda_{m})$ defines a path, i.e, the ordered product of non-linear and non-commuting operators

with $U[\emptyset]x=x$. A scattering transform along the path $p$ is defined as an integral, normalized by the response of a Dirac:

Each scattering coefficient $\overline{S}x(p)$ is invariant to a translation of $x$. We shall see that this transform has many similarities with the Fourier transform modulus, which is also translation invariant. However, a scattering is Lipschitz continuous to deformations as opposed to the Fourier transform modulus.

For classification, it is often better to compute localized descriptors which are invariant to translations smaller than a predefined scale $2^{J}$, while keeping the spatial variability at scales larger than $2^{J}$. This is obtained by localizing the scattering integral with a scaled spatial window $\phi_{2^{J}}(u)=2^{-2J}\phi(2^{-J}u)$. It defines a windowed scattering transform in the neighborhood of $u$:

with $S_{J}[\emptyset]x=x\star\phi_{2^{J}}$. For each path $p$, $S_{J}[p]x(u)$ is a function of the window position $u$, which can be subsampled at intervals proportional to the window size $2^{J}$. The averaging by $\phi_{2^{J}}$ implies that $S_{J}[p]x(u)$ is nearly invariant to translations $L_{c}x(u)=x(u-c)$ if $|c|\ll 2^{J}$. Section 3.1 proves that it is also stable relatively to deformations.

3 Scattering Convolution Network

If $p$ is a path of length $m$ then $S_{J}[p]x(u)$ is called scattering coefficient of order $m$ at the scale $2^{J}$. It is computed at the layer $m$ of a convolution network which is specified. For large scale invariants, several layers are necessary to avoid losing crucial information.

For appropriate wavelets, first order coefficients $S_{J}[\lambda_{1}]x$ are equivalent to SIFT coefficients . Indeed, SIFT computes the local sum of image gradient amplitudes among image gradients having nearly the same direction, in a histogram having $8$ different direction bins. The DAISY approximation shows that these coefficients are well approximated by $S_{J}[2^{j}r]x=|x\star\psi_{2^{j}r}|\star\phi_{2^{J}}(u)$ where $\psi_{2^{j}r}$ is the partial derivative of a Gaussian computed at the finest image scale $2^{j}$, for $8$ different rotations $r$. The averaging filter $\phi_{2^{J}}$ is a scaled Gaussian.

Partial derivative wavelets are well adapted to detect edges or sharp transitions but do not have enough frequency and directional resolution to discriminate complex directional structures. For texture analysis, many researchers have been using averaged wavelet coefficient amplitudes $|x\star\psi_{\lambda}|\star\phi_{J}(u)$, but calculated with a complex wavelet $\psi$ having a better frequency and directional resolution.

A scattering transform computes higher-order coefficients by further iterating on wavelet transforms and modulus operators. At a maximum scale $2^{J}$, wavelet coefficients are computed at frequencies $2^{j}\geq 2^{-J}$, and lower frequencies are filtered by $\phi_{2^{J}}(u)=2^{-2J}\phi(2^{-J}u)$. Since images are real-valued signals, it is sufficient to consider “positive” rotations $r\in G^{+}$ with angles in $[0,\pi)$:

with $\Lambda_{J}=\{\lambda=2^{j}r~{}:~{}r\in G^{+},j\geq-J\}$. For a Morlet wavelet $\psi$, the averaging filter $\phi$ is chosen to be a Gaussian. Let us emphasize that $2^{J}$ is a spatial scale variable whereas $\lambda=2^{j}r$ is assimilated to a frequency variable.

A wavelet modulus propagator keeps the low-frequency averaging and computes the modulus of complex wavelet coefficients:

Let $\Lambda_{J}^{m}$ be the set of all paths $p=(\lambda_{1},...,\lambda_{m})$ of length $m$. We denote $U[\Lambda_{J}^{m}]x=\{U[p]x\}_{p\in\Lambda_{J}^{m}}$ and $S_{J}[\Lambda_{J}^{m}]x=\{S_{J}[p]x\}_{p\in\Lambda_{J}^{m}}$. Since

and $S_{J}[p]x=U[p]x\star\phi_{2^{J}}$, it results that

This implies that $S_{J}[p]x$ can be computed along paths of length $m\leq m_{\max}$ by first calculating $U_{J}x=\{S_{J}[\emptyset]x\,,\,U[\Lambda_{J}^{1}]x\}$ and iteratively applying $U_{J}$ to each $U[\Lambda_{J}^{m}]x$ for increasing $m\leq m_{\max}$. This algorithm is illustrated in Figure 2.

A scattering transform thus appears to be a deep convolution network , with some particularities. As opposed to most convolution networks, a scattering network outputs coefficients $S_{J}[p]x$ at all layers $m\leq m_{\max}$, and not just at the last layer $m_{\max}$ . The next section proves that the energy of the deepest layer converges quickly to zero as $m_{\max}$ increases.

A second distinction is that filters are not learned from data but are predefined wavelets. Wavelets are stable with respect to deformations and provide sparse image representations. Stability to deformations is a strong condition which imposes a separation of the different image scales , hence the use of wavelets.

The modulus operator which recombines real and imaginary parts can be interpreted as a pooling function in the context of convolution networks. The averaging by $\phi_{2^{J}}$ at the output is also a pooling operator which aggregates coefficients to build an invariant. It has been argued that an average pooling loses information, which has motivated the use of other operators such as hierarchical maxima . The high frequencies lost by the averaging are recovered as wavelet coefficients in the next layers, which explains the importance of using a multilayer network structure. As a result, it only loses the phase of these wavelet coefficients. This phase may however be recovered from the modulus thanks to the wavelet transform redundancy. It has been proved that the wavelet-modulus operator $U_{J}x=\{x\star\phi_{2^{J}},\,|x\star\psi_{\lambda}|\}_{\lambda\in\Lambda_{J}}$ is invertible with a continuous inverse. It means that $x$ and hence the complex phase of each $x\star\psi_{\lambda}$ can be reconstructed. Although $U_{J}$ is invertible, the scattering transform is not exactly invertible because of instabilities. Indeed, applying $U_{J}$ in (7) for $m\leq m_{\max}$ computes all $S_{J}[\Lambda_{J}^{m}]x$ for $m\leq m_{\max}$ but also the last layer of internal network coefficients $U[\Lambda_{J}^{m_{\max}+1}]x$. The next section proves that $U[\Lambda_{J}^{m_{\max}+1}]x$ can be neglected because its energy converges to zero as $m_{\max}$ increases. However, this introduces a small error which accumulates when iterating on $U_{J}^{-1}$.

Figure 4 shows the Fourier transform of two images, and the amplitude of their scattering coefficients of orders $m=1$ and $m=2$, at a maximum scale $2^{J}$ equal to the image size. A scattering coefficient over a quadrant $\Omega[2^{j_{1}}r_{1}]$ gives an approximation of the Fourier transform energy over the support of $\hat{\psi}_{2^{j_{1}}r_{1}}$. Although the top and bottom images are very different, they have same order $m=1$ scattering coefficients. Here, first-order coefficients are not sufficient to discriminate between two very different images. However, coefficients of order $m=2$ succeed in discriminating between the two images. The top image has wavelet coefficients which are much more sparse than the bottom image. As a result, Section 3.1 shows that second-order scattering coefficients have a larger amplitude. Higher-order coefficients are not displayed because they have a negligible energy as explained in Section 3.

Scattering Properties

A convolution network is highly non-linear, which makes it difficult to understand how the coefficient values relate to the signal properties. For a scattering network, Section 3.1 analyzes the coefficient properties and optimizes the network architecture. For texture analysis, the scattering transform of stationary processes is studied in Section 3.2. The regularity of scattering coefficients can be exploited to reduce the size of a scattering representation, by using a cosine transform, as shown in Section 3.3. Finally, Section 3.4 provides a fast computational algorithm.

A windowed scattering $S_{J}$ is computed with a cascade of wavelet modulus operators $U_{J}$, and its properties thus depend upon the wavelet transform properties. Conditions are given on wavelets to define a scattering transform which is contracting and preserves the signal norm. This analysis shows that $\|S_{J}[p]x\|$ decreases quickly as the length of $p$ increases, and is non-negligible only over a particular subset of frequency-decreasing paths. Reducing computations to these paths defines a convolution network with much fewer internal and output coefficients.

then applying the Plancherel formula proves that $W_{J}x=\{x\star\phi_{J}\,,\,x\star\psi_{\lambda}\}_{\lambda\in\Lambda_{J}}$ satisfies

with $\|W_{J}x\|^{2}=\|x\star\phi_{J}\|^{2}+\sum_{\lambda\in\Lambda_{J}}\|x\star\psi_{\lambda}\|^{2}$. In the following we suppose that $\epsilon<1$ and hence that the wavelet transform is a contracting and invertible operator, with a stable inverse. If $\epsilon=0$ then $W_{J}$ is unitary. The Morlet wavelet $\psi$ in Figure 1 satisfies (8) with $\epsilon=0.25$, together with $\phi(u)=C\exp(-|u|^{2}/(2\sigma_{0}^{2}))$ with $\sigma_{0}=0.7$ and $C$ adjusted so that $\int\phi(u)\,du=1$. These functions are used in all classification applications. Rotated and dilated cubic spline wavelets are constructed in to satisfy (8) with $\epsilon=0$.

The modulus is contracting in the sense that $||a|-|b||\leq|a-b|$. Since $U_{J}=\{x\star\phi_{J}\,,\,|x\star\psi_{\lambda}|\}_{\lambda\in\Lambda_{J}}$ is obtained with a wavelet transform $W_{J}$ followed by modulus, which are both contractive, it is also contractive:

If $W_{J}$ is unitary then $U_{J}$ also preserves the signal norm $\|U_{J}x\|=\|x\|$.

If $W_{J}$ is unitary, $\epsilon=0$ in (9) and for appropriate wavelets, it is proved in that

This result uses the fact that $U_{J}$ preserves the signal norm and that $U_{J}\,U[\Lambda_{J}^{m}]x=\{S_{J}[\Lambda_{J}^{m}]x\,,\,U[\Lambda_{J}^{m+1}]x\}$. Proving (10) is thus equivalent to prove that the energy of the last network layer converges to zero when $m_{\max}$ increases

This result is also important for numerical applications because it explains why the network depth can be limited with a negligible loss of signal energy.

The scattering energy conservation also provides a relation between the network energy distribution and the wavelet transform sparsity. For $p=(\lambda_{1},...,\lambda_{m})$, we denote $p+\lambda=(\lambda,\lambda_{1},...,\lambda_{m})$. Applying (10) to $U[\lambda]x=|x\star\psi_{\lambda}|$ instead of $x$, and separating the first term for $m=0$ yields

The energy conservation (10) is proved by showing that the scattering energy $\|U[p]x\|^{2}$ propagates towards lower frequencies as the length of $p$ increases. This energy is thus ultimately captured by the low-pass filter $\phi_{2^{J}}$ which outputs $S_{J}[p]x=U[p]x\star\phi_{2^{J}}$. This property requires that $x\star\psi_{\lambda}$ has a lower-frequency envelope $|x\star\psi_{\lambda}|$. It is valid if $\psi(u)=e^{i\eta.u}\,\theta(u)$ where $\theta$ is a low-pass filter. To verify this property, we write $x\star\psi_{\lambda}(u)=e^{i\lambda\xi.u}\,x_{\lambda}(u)$ with

This signal is filtered by the dilated and rotated low-pass filter $\theta_{\lambda}$ whose Fourier transform is $\hat{\theta}_{\lambda}(\omega)=\theta(\lambda^{-1}\omega)$. So $|x\star\psi_{\lambda}(u)|=|x_{\lambda}(u)|$ is the modulus of a regular function and is therefore mostly regular. This result is not valid if $\psi$ is a real because $|x\star\psi_{\lambda}|$ is singular at each zero-crossing of $x\star\psi_{\lambda}(u)$.

The modulus appears as a non-linear “demodulator” which projects wavelet coefficients to lower frequencies. If $\lambda=2^{j}r$ then $|x\star\psi_{\lambda}(u)|\star\psi_{\lambda^{\prime}}$ for $\lambda^{\prime}=2^{j^{\prime}}r^{\prime}$ is non-negligible only if $\psi_{\lambda^{\prime}}$ is located at low frequencies and hence if $2^{j^{\prime}}<2^{j}$. Iterating on wavelet modulus operators thus propagates the scattering energy along frequency-decreasing paths $p=(2^{j_{1}}r_{1},...,2^{j_{m}}r_{m})$ where $2^{j_{k}}\leq 2^{j_{k-1}}$, for $1\leq k<m$. Scattering coefficients along other paths have a negligible energy. Over the CalTech101 images database, Table I shows that over $99\%$ of the scattering energy is concentrated along frequency-decreasing paths of length $m\leq 3$. Numerically, it is therefore sufficient to compute the scattering transform along this subset of frequency-decreasing paths. It defines a much smaller convolution network. Section 3.4 shows that the resulting coefficients are computed with $O(N\log N)$ operations.

For classification applications, one of the most important properties of a scattering transform is its stability to deformations $L_{\tau}x(u)=x(u-\tau(u))$, because wavelets are stable to deformations and the modulus commutes with $L_{\tau}$. Let $\|\tau\|_{\infty}=\sup_{u}|\tau(u)|$ and $\|\nabla\tau\|_{\infty}=\sup_{u}|\nabla\tau(u)|<1$. If $S_{J}$ is computed on paths of length $m\leq m_{\max}$ then it is proved in that for signals $x$ of compact support

with a second order Hessian term which is negligible if $\tau(u)$ is regular. If $2^{J}\geq\|\tau\|_{\infty}/\|\nabla\tau\|_{\infty}$ then the translation term can be neglected and the transform is Lipschitz continuous to deformations:

2 Scattering Stationary Processes

Image textures can be modeled as realizations of stationary processes $X(u)$. We denote the expected value of $X$ by $E(X)$, which does not depend upon $u$. The Fourier spectrum $\widehat{R}X(\omega)$ is the Fourier transform of the autocorrelation

Despite the importance of spectral methods, the Fourier spectrum is often not sufficient to discriminate image textures because it does not take into account higher-order moments. Figure 5 shows very different textures having same second-order moments. A scattering representation of stationary processes includes second order and higher-order moment descriptors of stationary processes, which discriminates between such textures.

If $X(u)$ is stationary then $U[p]X(u)$ remains stationary because it is computed with a cascade of convolutions and modulus operators which preserve stationarity. Its expected value thus does not depend upon $u$ and defines the expected scattering transform:

A windowed scattering gives an estimator of $\overline{S}X(p)$, calculated from a single realization of $X$, by averaging $U[p]X$ with $\phi_{2^{J}}$:

Since $\int\phi_{2^{J}}(u)\,du=1$, this estimator is unbiased: $E(S_{J}[p]X)=E(U[p]X)$.

For appropriate wavelets, it is also proved that

Replacing $X$ by $X\star\psi_{\lambda}$ implies that

These expected squared wavelet coefficients can also be written as a filtered integration of the Fourier power spectrum $\widehat{R}X(\omega)$

These two equations prove that summing scattering coefficients recovers the power spectrum integral over each wavelet frequency support, which only depends upon second-order moments. However, one can also show that scattering coefficients $\overline{S}X(p)$ depend upon moments of $X$ up to the order $2^{m}$ if $p$ has a length $m$. Scattering coefficients can thus discriminate textures having same second-order moments but different higher-order moments. This is illustrated using the two textures in Figure 5, which have the same power spectrum and hence same second order moments. Scattering coefficients $S_{J}[p]X$ are shown for $m=1$ and $m=2$ with the frequency tiling illustrated in Figure 3. The ability to discriminate the top process $X_{1}$ from the bottom process $X_{2}$ is measured by a scattering distance normalized by the variance:

For $m=1$, scattering coefficients mostly depend upon second-order moments and are thus nearly equal for both textures. One can indeed verify numerically that $\rho(1)=1$ so both textures can not be distinguished using first order scattering coefficients. On the contrary, scattering coefficients of order $2$ are highly dissimilar because they depend on moments up to order $4$, and $\rho(2)=5$.

For a large class of ergodic processes including most image textures, it is observed numerically that the total scattering variance $\sum_{p\in{\cal P}_{J}}E(|S_{J}[p]X-\overline{S}X(p)|^{2})$ decreases to zero when $2^{J}$ increases. Table II shows the decay of the total scattering variance, computed on average over the Brodatz texture dataset. Since $E(|S_{J}[p]X|^{2})=E(S_{J}[p]X)^{2}+E(|S_{J}[p]X-E(S_{J}[p]X)|^{2})$ and $E(S_{J}[p]X)=\overline{S}X(p)$, it results from the energy conservation (15) that the expected scattering transform also satisfies

Table III gives the percentage of expected scattering energy $\sum_{p\in\Lambda_{\infty}^{m}}|\overline{S}X(p)|^{2}$ carried by paths of length $m$, for textures in the Brodatz database. Most of the energy is concentrated in paths of length $m\leq 3$.

3 Cosine Scattering Transform

Natural images have scattering coefficients $S_{J}[p]X(u)$ which are correlated across paths $p=(2^{j_{1}}r_{1},...,2^{j_{m}}r_{m})$, at any given position $u$. The strongest correlation is between paths of same length. For each $m$, scattering coefficients are decorrelated in a Karhunen-Loève basis which diagonalizes their covariance matrix. Figure 6 compares the decay of the sorted variances $E(|S_{J}[p]X-E(S_{J}[p]X)|^{2})$ and the variance decay in the Karhunen-Loève basis computed on paths of length $m=1$, and on paths of length $m=2$, over the Caltech image dataset with a Morlet wavelet. The variance decay is much faster in the Karhunen-Loève basis, which shows that there is a strong correlation between scattering coefficients of same path length.

A change of variables proves that a rotation and scaling $X_{2^{l}r}(u)=X(2^{-l}ru)$ produces a rotation and inverse scaling on the path variable $p=(2^{j_{1}}r_{1},...,2^{j_{m}}r_{m})$:

If images are randomly rotated and scaled by $2^{l}r^{-1}$ then the path $p$ is randomly rotated and scaled . In this case, the scattering transform has stationary variations along the scale and rotation variables. This suggests approximating the Karhunen-Loève basis by a cosine basis along these variables. Let us parameterize each rotation $r$ by its angle $\theta\in[0,2\pi]$. A path $p$ is then parameterized by $([j_{1},\theta_{1}],...,[j_{m},\theta_{m}])$.

4 Fast Scattering Computations

Section 3.1 shows that the scattering energy is concentrated along frequency-decreasing paths $p=(2^{j_{k}}r_{k})_{k}$ satisfying $2^{-J}\leq 2^{j_{k+1}}<2^{j_{k}}$. If the wavelet transform is computed along $C$ directions then the total number of frequency-decreasing paths of length $m$ is $C^{m}\binom{J}{m}$. Since $\phi_{2^{J}}$ is a low-pass filter, $S_{J}[p]x(u)=U[p]x\star\phi_{2^{J}}(u)$ can be uniformly sampled at intervals $\alpha 2^{J}$, with $\alpha=1$ or $\alpha=1/2$. If $x(n)$ is a discrete image with $N$ pixels, then each $S_{J}[p]x$ has $2^{-2J}\alpha^{-2}N$ coefficients. The scattering representation along all frequency-decreasing paths of length at most $m$ thus has a total number of coefficients equal to $N_{J}=N\alpha^{-2}2^{-2J}\sum_{q=0}^{m}C^{q}\binom{J}{q}$. This reduced scattering representation is computed by a cascade of convolutions, modulus, and sub-samplings, with $O(N\log N)$ operations. The final DCT transform further compresses the resulting representation.

Let us recall from Section 2.3 that scattering coefficients are computed by iteratively applying the one-step propagator $U_{J}$. To compute subsampled scattering coefficients along frequency-decreasing paths, this propagator is truncated. For any scale $2^{k}$, ${U}_{k,J}$ transforms a signal $x(2^{k}\alpha n)$ into

The algorithm computes subsampled scattering coefficients by iterating on this propagator.

If $x$ is a signal of size $P$ then FFT’s compute $U_{k,J}x$ with $O(P\log P)$ operations. A reduced scattering transform thus computes its $N_{J}=N\alpha^{-2}2^{-2J}\sum_{m=0}^{m_{\max}}C^{m}\binom{J}{m}$ coefficients with $O(N_{J}\log N)$ operations. If $m_{\max}=2$ then $N_{J}=N\alpha^{-2}2^{-2J}(CJ+C^{2}J(J-1)/2)$. It decreases exponentially when the scale $2^{J}$ increases.

Numerical computations in this paper are performed by rotating wavelets along $C=6$ directions, for scattering representations of maximum order $m_{\max}=2$. The resulting size of a reduced cosine scattering representation has at most three times as many coefficients as a dense SIFT representation. SIFT represents small blocks of $4^{2}$ pixels with $8$ coefficients. A cosine scattering representation represents each image block of $2^{2J}$ pixels by $N_{J}2^{2J}/(2N)=(CJ+C^{2}J(J-1)/2)/2$ coefficients, which is equal to $24$ for $C=6$ and $J=2$. The cosine scattering transform is thus three times the size of SIFT for $J=2$, but as $J$ increases, the relative size decreases. If $J=3$ then the size of a cosine scattering representation is twice the size of a SIFT representation but for $J=7$ it is about $20$ times smaller.

Classification Using Scattering Vectors

A scattering transform eliminates the image variability due to translation and is stable to deformations. The resulting classification properties are studied with a PCA and an SVM classifier applied to scattering representations computed with a Morlet wavelet. State-of-the-art results are obtained for hand-written digit recognition and for texture discrimination.

Although discriminant classifiers such as SVM have better asymptotic properties than generative classifiers , the situation can be inverted for small training sets. We introduce a simple robust generative classifier based on affine space models computed with a PCA. Applying a DCT on scattering coefficients has no effect on any linear classifier because it is a linear orthogonal transform. However, keeping the $50$% lower frequency cosine scattering coefficients reduces computations and has a negligible effect on classification results. The classification algorithm is described directly on scattering coefficients to simplify explanations. Each signal class is represented by a random vector $X_{k}$, whose realizations are images of $N$ pixels in the class.

Let $E(S_{J}X)=\{E({S}_{J}[p]X(u))\}_{p,u}$ be the family of $N_{J}$ expected scattering values, computed along all frequency-decreasing paths of length $m\leq m_{\max}$ and all subsampled positions $u=\alpha 2^{J}n$. The difference ${S}_{J}X_{k}-E({S}_{J}X_{k})$ is approximated by its projection in a linear space of low dimension $d\ll N_{J}$. The covariance matrix of ${S}_{J}X_{k}$ is a matrix of size $N_{J}^{2}$. Let ${\mathbf{V}_{d,k}}$ be the linear space generated by the $d$ PCA eigenvectors of this covariance matrix having the largest eigenvalues. Among all linear spaces of dimension $d$, this is the space which approximates ${S}_{J}X_{k}-E({S}_{J}X_{k})$ with the smallest expected quadratic error. This is equivalent to approximating ${S}_{J}X_{k}$ by its projection on an affine approximation space:

The resulting classifier associates a signal $X$ to the class $\hat{k}$ which yields the best approximation space:

The minimization of this distance has similarities with the minimization of a tangential distance in the sense that we remove the principal scattering directions of variabilities to evaluate the distance. However it is much simpler since it does not evaluate a tangential space which depends upon ${S}_{J}x$. Let $\mathbf{V}_{d,k}^{\perp}$ be the orthogonal complement of $\mathbf{V}_{d,k}$ corresponding to directions of lower variability. This distance is also equal to the norm of the difference between ${S}_{J}x$ and the average class “template” $E({S}_{J}X_{k})$, projected in $\mathbf{V}_{d,k}^{\perp}$:

Minimizing the affine space approximation error is thus equivalent to finding the class centroid $E({S}_{J}X_{k})$ which is the closest to ${S}_{J}x$, without taking into account the first $d$ principal variability directions. The $d$ principal directions of the space $\mathbf{V}_{d,k}$ result from deformations and from structural variability. The projection $P_{\mathbf{A}_{d,k}}({S}_{J}x)$ is the optimum linear prediction of ${S}_{J}x$ from these $d$ principal modes. The selected class has the smallest prediction error.

This affine space selection is effective if ${S}_{J}X_{k}-E({S}_{J}X_{k})$ is well approximated by a projection in a low-dimensional space. This is the case if realizations of $X_{k}$ are translations and limited deformations of a single template. Indeed, the Lipschitz continuity condition implies that small deformations are linearized by the scattering transform. Hand-written digit recognition is an example. This is also valid for stationary textures where ${S}_{J}X_{k}$ has a small variance, which can be interpreted as structural variability.

The dimension $d$ must be adjusted so that ${S}_{J}X_{k}$ has a better approximation in the affine space ${\bf A}_{d,k}$ than in affine spaces ${\bf A}_{d,k^{\prime}}$ of other classes $k^{\prime}\neq k$. This is a model selection problem, which requires to optimize the dimension $d$ in order to avoid over-fitting .

The invariance scale $2^{J}$ must also be optimized. When the scale $2^{J}$ increases, translation invariance increases but it comes with a partial loss of information which brings the representations of different signals closer. One can prove that for any $x$ and $x^{\prime}$

When $2^{J}$ goes to infinity, this scattering distance converges to a non-zero value. To classify deformed templates such as hand-written digits, the optimal $2^{J}$ is of the order of the maximum pixel displacements due to translations and deformations. In a stochastic framework where $x$ and $x^{\prime}$ are realizations of stationary processes, $S_{J}x$ and $S_{J}x^{\prime}$ converge to the expected scattering transforms $\overline{S}x$ and $\overline{S}x^{\prime}$. In order to classify stationary processes such as textures, the optimal scale is the maximum scale equal to the image width, because it minimizes the variance of the windowed scattering estimator.

A cross-validation procedure is used to find the dimension $d$ and the scale $2^{J}$ which yield the smallest classification error. This error is computed on a subset of the training images, which is not used to estimate the covariance matrix for the PCA calculations.

As in the case of SVM, the performance of the affine PCA classifier can be improved by equalizing the descriptor space. Table I shows that scattering vectors have unequal energy distribution along its path variables, in particular as the order varies. A robust equalization is obtained by re-normalizing each $S_{J}[p]X(u)$ by the maximum \|S_{J}[p]X_{i}\|=\Big{(}\sum_{u}|S_{J}[p]X_{i}(u)|^{2}\Big{)}^{1/2} over all training signals $X_{i}$:

To simplify notations, we still write $S_{J}X$ for this normalized scattering vector.

Affine space scattering models can be interpreted as generative models computed independently for each class. As opposed to discriminative classifiers such as SVM, they do not estimate cross-terms between classes, besides from the choice of the model dimensionality $d$. Such estimators are particularly effective for small number of training samples per class. Indeed, if there are few training samples per class then variance terms dominate bias errors when estimating off-diagonal covariance coefficients between classes .

An affine space approximation classifier can also be interpreted as a robust quadratic discriminant classifier obtained by coarsely quantizing the eigenvalues of the inverse covariance matrix. For each class, the eigenvalues of the inverse covariance are set to in ${\bf V}_{d,k}$ and to $1$ in ${\bf V}^{\perp}_{d,k}$, where $d$ is adjusted by cross-validation. This coarse quantization is justified by the poor estimation of covariance eigenvalues from few training samples. These affine space models will typically be applied to distributions of scattering vectors having non-Gaussian distributions, where a Gaussian Fisher discriminant can lead to important errors.

2 Handwritten Digit Recognition

The MNIST database of hand-written digits is an example of structured pattern classification, where most of the intra-class variability is due to local translations and deformations. It comprises at most $60000$ training samples and $10000$ test samples. If the training dataset is not augmented with deformations, the state of the art was achieved by deep-learning convolutional networks , deformation models , and dictionary learning . These results are improved by a scattering classifier.

All computations are performed on the reduced cosine scattering representation described in Section 3.3, which keeps the lower-frequency half of the coefficients. Table IV computes classification errors on a fixed set of test images, depending upon the size of the training set, for different representations and classifiers. The affine space selection of section 4.1 is compared with an SVM classifier using RBF kernels, which are computed using Libsvm , and whose variance is adjusted using standard cross-validation over a subset of the training set. The SVM classifier is trained with a renormalization which maps all coefficients to $$. The PCA classifier is trained with the renormalisation (19). The first two columns of Table IV show that classification errors are much smaller with an SVM than with the PCA algorithm if applied directly on the image. The 3rd and 4th columns give the classification error obtained with a PCA or an SVM classification applied to the modulus of a windowed Fourier transform. The spatial size$2^{J}$of the window is optimized with a cross-validation which yields a minimum error for$2^{J}=8$. It corresponds to the largest pixel displacements due to translations or deformations in each class. Removing the complex phase of the windowed Fourier transform yields a locally invariant representation but whose high frequencies are unstable to deformations, as explained in Section 2.1. Suppressing this local translation variability improves the classification rate by a factor$3$for a PCA and by almost$2$ for an SVM. The comparison between PCA and SVM confirms the fact that generative classifiers can outperform discriminative classifiers when training samples are scarce . As the training set size increases, the bias-variance trade-off turns in favor of the richer SVM classifiers, independently of the descriptor.

Columns 6 and 8 give the PCA classification result applied to a windowed scattering representation for $m_{\max}=1$ and $m_{\max}=2$. The cross validation also chooses $2^{J}=8$. For the digit ‘3’, Figure 7 displays the 4-by-4 array of normalized scattering vectors. For each $u=2^{J}(n_{1},n_{2})$ with $1\leq n_{i}\leq 4$, the scattering vector $S_{J}[p]X(u)$ is displayed for paths of length $m=1$ and $m=2$, as circular frequency energy distributions following Section 2.3.

Increasing the scattering order from $m_{\max}=1$ to $m_{\max}=2$ reduces errors by about $30$%, which shows that second order coefficients carry important information even at a relatively small scale $2^{J}=8$. However, third order coefficients have a negligible energy and including them brings marginal classification improvements, while increasing computations by an important factor. As the learning set increases in size, the classification improvement of a scattering transform increases relatively to windowed Fourier transform because the classification is able to incorporate more high frequency structures, which have deformation instabilities in the Fourier domain as opposed to the scattering domain.

Table IV also shows that below $5\cdot\,10^{3}$ training samples, the scattering PCA classifier improves results of a deep-learning convolutional networks, which learns all filter coefficients with a back-propagation algorithm . As more training samples are available, the flexibility of the SVM classifier brings an improvement over the more rigid affine classifier, yielding a $0.43\%$ error rate on the original dataset, thus improving upon previous state of the art methods.

To evaluate the precision of the affine space model, we compute the relative affine approximation error, averaged over all classes:

For any ${S}_{J}X_{k}$, we also calculate the minimum approximation error produced by another affine model $A_{d,k^{\prime}}$ with $k^{\prime}\neq k$:

For a scattering representation with $m_{\max}=2$, Table V gives the dimension $d$ of affine approximation spaces optimized with a cross validation, with the corresponding values of $\sigma_{d}^{2}$ and $\lambda_{d}$. When the training set size increases, the model dimension $d$ increases because there are more samples to estimate each intra-class covariance matrix. The approximation model becomes more precise so $\sigma_{d}^{2}$ decreases and the relative approximation error $\lambda_{d}$ produced by wrong classes increases. This explains the reduction of the classification error rate observed in Table IV as the training size increases.

The US-Postal Service is another handwritten digit dataset, with 7291 training samples and 2007 test images $16\times 16$ pixels. The state of the art is obtained with tangent distance kernels . Table VI gives results obtained with a scattering transform with the PCA classifier for $m_{\max}=1,2$. The cross-validation sets the scattering scale to $2^{J}=8$. As in the MNIST case, the error is reduced when going from $m_{\max}=1$ to $m_{\max}=2$ but remains stable for $m_{\max}=3$. Different renormalization strategies can bring marginal improvements on this dataset. If the renormalization is performed by equalizing using the standard deviation of each component, the classification error is ${2.3\%}$ whereas it is $2.6\%$ if the supremum is normalized.

The scattering transform is stable but not invariant to rotations. Stability to rotations is demonstrated over the MNIST database in the setting defined in . A database with 12000 training samples and 50000 test images is constructed with random rotations of MNIST digits. The PCA affine space selection takes into account the rotation variability by increasing the dimension $d$ of the affine approximation space. This is equivalent to projecting the distance to the class centroid on a smaller orthogonal space, by removing more principal components. The error rate in Table VII is much smaller with a scattering PCA than with a convolution network . Much better results are obtained for a scattering with $m_{\max}=2$ than with $m_{\max}=1$ because second order coefficients maintain enough discriminability despite the removal of a larger number $d$ of principal directions. In this case, $m_{\max}=3$ marginally reduces the error.

Scaling invariance is studied by introducing a random scaling factor uniformly distributed between $1/\sqrt{2}$ and $\sqrt{2}$. In this case, the digit ‘9’ is removed from the database as to avoid any indetermination with the digit ‘6’ when rotated. The training set has $9000$ samples ($1000$ samples per class). Table VIII gives the error rate on the original MNIST database and including either rotation, scaling, or both in the training and testing samples. Scaling has a smaller impact on the error rate than rotating digits because scaled scattering vectors span an invariant linear space of lower dimension. Second-order scattering outperforms first-order scattering, and the difference becomes more significant when rotation and scaling are combined, because it provides interaction coefficients which are discriminative even in presence of scaling and rotation variability.

3 Texture Discrimination

Visual texture discrimination remains an outstanding image processing problem because textures are realizations of non-Gaussian stationary processes, which cannot be discriminated using the power spectrum. Depending on the imaging conditions, textures undergo transformations due to illumination, rotation, scaling or more complex deformations when mapped on three-dimensional surfaces. The affine PCA space classifier removes most of the variability of $S_{J}X-E\{S_{J}X\}$ within each class. This variability is due to the residual stochastic variability which decays as $J$ increases and to variability due to illumination, rotation and perspective effects.

Texture classification is tested on the CUReT texture database , which includes 61 classes of image textures of $N=200^{2}$ pixels. Each texture class gives images of the same material with different pose and illumination conditions. Specularities, shadowing and surface normal variations make classification challenging. Pose variation requires global rotation and illumination invariance. Figure 8 illustrates the large intra-class variability, after a normalization of the mean and variance of each textured image.

Table IX compares error rates obtained with different classifiers. The database is randomly split into a training and a testing set, with 46 training images for each class as in . Results are averaged over 10 different splits. A PCA affine space classifier applied directly on the image yields a large classification error of $17\%$. To estimate the Fourier spectrum, windowed Fourier transforms are computed over half-overlapping windows of size $2^{J}$, and their squared modulus is averaged over the whole image. This averaging is necessary to reduce the spectrum estimator variance, which does not decrease when the window size $2^{J}$ increases. The cross-validation sets the optimal window scale to $2^{J}=32$, whereas images have a width of $200$ pixels. The error drops to 1%. This simple Fourier spectrum yields a smaller error than previously reported state-of-the-art methods. SVM’s applied to a dictionary of textons yield an error rate of 1.53% , whereas an optimized Markov Random Field model computed with image patches achieves an error of 2.46%.

For the scattering PCA classifier, the cross validation chooses an optimal scale $2^{J}$ equal to the image width to reduce the scattering estimation variance. Indeed, contrarly to a power spectrum estimation, the variance of the scattering vector decreases when $2^{J}$ increases. Figure 9 displays the scattering coefficients $S_{J}[p]X$ of order $m=1$ and $m=2$ of a CureT textured image $X$. When $m_{max}=1$, the error drops to 0.5%, although first-order scattering coefficients essentially depend upon second order moments as the Fourier spectrum. This is probably due to the fact that image textures have a spectrum which typically decays like $|\omega|^{-\alpha}$. For such spectrum, an estimation over dyadic frequency bands provide a better bias versus variance trade-off than a windowed Fourier spectrum . For $m_{max}=2$, the error further drops to 0.2%. Indeed, scattering coefficients of order $m=2$ depend upon moments of order $4$, which are necessary to differentiate textures having same second order moments as in Figure 5. The dimension of the affine approximation space model is $d=16$, the intra-class normalized approximation error is $\sigma^{2}_{d}=2.5\cdot 10^{-1}$ and the separation ratio is $\lambda_{d}=3$ for $m_{\max}=2$.

The PCA classifier provides a partial rotation invariance by removing principal components. It averages scattering coefficients along path rotation parameters, which comes with a loss of discriminability. However, a more efficient rotation invariant texture classification is obtained by cascading this translation invariant scattering with a second rotation invariant scattering . It transforms each layer of the translation invariant scattering network with new wavelet convolutions along rotation parameters, followed by modulus and average pooling operators, which are cascaded. A combined translation and rotation scattering yields a translation and rotation invariant representation which is stable to deformations .

Conclusion

A wavelet scattering transform computes a translation invariant representation, which is stable to deformation, using a deep convolution network architecture. The first layer outputs SIFT-type descriptors, which are not sufficiently informative for large-scale invariance. Classification performance is improved by adding other layers providing complementary information. A reduced cosine scattering transform is at most three times larger than a SIFT descriptor and computed with $O(N\log N)$ operations.

State-of-the-art classification results are obtained for handwritten digit recognition and texture discrimination, with an SVM or a PCA classifier. If the data set has other sources of variability due to the action of other finite Lie groups such as rotations, then this variability can be eliminated with an invariant scattering computed by cascading wavelet transforms defined on these groups .

However, signal classes may also include complex sources of variability that can not be approximated by the action of a finite group, as in CalTech101 or Pascal databases. This variability must be taken into account by unsupervised optimizations of the representations from the training data. Deep convolution networks which learn filters from the data have the flexibility to adapt to such variability, but learning translation invariant filters is not necessary. A wavelet scattering transform can be used on the first two network layers, while learning the next layer filters applied to scattering coefficients. Similarly, bag-of-features unsupervised algorithms applied to SIFT can potentially be improved upon by replacing SIFT descriptors by wavelet scattering vectors.

References