Classification with Scattering Operators

Introduction

Locally invariant image descriptors such as SIFT provide efficient image representations for image classification and registration . These feature vectors as well as multiscale texture descriptors can be computed with a spatial averaging of wavelet coefficient amplitudes. The averaging reduces the feature variability and provides local translation invariance, but it also reduces information.

Scattering operators recover the lost high frequencies and retransform them into co-occurrence coefficients at multiple scales and orientations. They provide much richer descriptors of complex structures such as corners, junctions and multiscale texture variations. These coefficients are locally translation invariant and they linearize small deformations. They are computed with a convolution network which cascades contractive wavelet transforms and modulus operators . Scattering operators provide new representations of stationary image textures, which can discriminate texture having the same power spectrum.

The scattering transform of a class of signals is approximated by an affine space computed with a PCA. Images are classified by selecting a best approximation space model for their scattering transform. State of the art results are obtained for hand-written digit recognition and for texture discrimination, with important rotation and illumination variability, and small training sets.

Section 2.1 reviews the relations between wavelet transforms and computer vision descriptors. Section 2.2 introduces scattering image representations. Classification by scattering model selection is introduced in Section 3, with numerical results. Softwares are available at www.cmap.polytechnique.fr/scattering.

Scattering

A scattering transform computes local image descriptors with a cascade of wavelet decompositions, complex modulus and a local averaging. The resulting scattering representation is locally invariant to translations. It includes coefficients which are similar to SIFT descriptors, together with co-occurrences coefficients at multiple scales and orientations.

Image feature vectors such as SIFT and multiscale Gabor textons are obtained by averaging the amplitude of wavelet coefficients, calculated with directional wavelets. Writing these feature vectors as wavelet coefficients helps to understand and to improve their properties.

The directional wavelet transform of $f$ at a position $x$ for scales $2^{j}<2^{J}$ is a vector of coefficients

where $\phi_{J}(x)=2^{-2J}\phi(2^{-J}x)$ is a low-pass filter which carries the low frequencies of $f$ above the scale $2^{J}$: $\int\phi(x)dx=1$. Let $|W_{J}f(x)|^{2}$ be the Euclidean norm of this vector which sums the square of its coordinates. Let $\hat{f}(\omega)$ be the Fourier transform of $f$. If wavelets satisfy

and this inequality is an equality if (2) is an equality. The wavelet transform is then contractive and potentially unitary.

Many standard image feature vectors are obtained by averaging wavelet coefficient amplitudes. SIFT coefficients are obtained from histograms of image gradients calculated at a fine scale $2^{j}$. A histogram bin indexed by $\gamma\in\Gamma$ stores the local sum of the amplitudes of all gradient vectors whose orientations are close to $\gamma$. Several authors observed that approximate SIFT feature vectors are computed more efficiently by averaging directly the partial derivative amplitudes of $f$ along the $K$ directions $\gamma\in\Gamma$, with a low-pass filter $\phi_{J}$. These averaged partial derivative amplitudes can be written as averaged wavelet coefficients

with a partial derivative wavelet $\psi(x)={\partial g(x)}/{\partial x_{1}}$, with $g(x)=e^{-|x|^{2}/2}$ and $x=(x_{1},x_{2})$. These averaged wavelet coefficients are nearly invariant to translations or deformations which are small relatively to $2^{J}$.

Partial derivative wavelets are well adapted to detect edge type elements, but these wavelets do not have enough frequency and directional resolution to discriminate more complex structures appearing in textures. For texture analysis, wavelets with a better frequency localization are often used . Complex Gabor functions are examples of such directional wavelets obtained by modulating a Gaussian window at a frequency $\xi$:

For stationary textures, $|f\star\psi_{j,\gamma}|\star\phi_{J}(x)$ has a reduced stochastic variability because of the averaging kernel $\phi_{J}$.

2 Scattering Coefficients

The local translation invariance and variability reduction of SIFT descriptors and multiscale textons is obtained by averaging. Scattering operators restore part of the information lost by this averaging with co-occurrence coefficients having similar invariance properties.

The wavelet transform (1) shows that high frequencies eliminated in $|f\star\psi_{j_{1},\gamma_{1}}|\star\phi_{J}$ by the convolution with $\phi_{J}$ are recovered by convolutions with wavelets $|f\star\psi_{j_{1},\gamma_{1}}|\star\psi_{j_{2},\gamma_{2}}$ at scales $2^{j_{2}}<2^{J}$. To become insensitive to local translation and reduce the variability of these coefficients, their complex phase is removed by a modulus, and it is averaged by $\phi_{J}$:

These are called scattering coefficients because they result from all interferences of $f$ with two wavelets . They give co-occurrence information in $f$ for any pair of scales $2^{j_{1}}$, $2^{j_{2}}$ and any two directions ${\gamma_{1}}$ and $\gamma_{2}$. This can distinguish corners and junctions from edges and it characterizes texture structures. Coefficients are only calculated for $2^{j_{2}}<2^{j_{1}}$ because one can show that $|f\star\psi_{j_{1},\gamma_{1}}|\star\psi_{j_{2},\gamma_{2}}$ is negligible at scales $2^{j_{2}}\geq 2^{j_{1}}$.

The convolution with $\phi_{J}$ removes high frequencies and thus yields second order coefficients that are locally translation invariant. High frequencies can again be restored by finer scale wavelet coefficients, which are regularized by averaging their amplitude with $\phi_{J}$. Applying iteratively this procedure $q$ times yields a vector of coefficients at each $x$:

This vector has $K^{q}\binom{J}{q}$ scattering coefficients, computing interactions between $f$ and the successive wavelets $\psi_{j_{1},\gamma_{1}}\,$…$\,\psi_{j_{q},\gamma_{q}}$. A scattering vector aggregates all these coefficients up to a maximum order $q\leq m$:

and the first coefficient is the signal average $S_{0,J}f(x)=f\star\phi_{J}(x)$. The scattering vector size is $\sum_{q=0}^{m}K^{q}\binom{J}{q}$. After convolution with $\phi_{J}$ the output is subsampled at intervals $2^{J}$. If $f(n)$ is an image of $N$ pixels, this uniform sampling yields a scattering representation $S_{J}f(2^{J}n)$ including a total of $N_{J}=2^{-2J}N\sum_{q=0}^{m}K^{q}\binom{J}{q}$ coefficients.

A scattering vector is computed with a cascade of convolutions and modulus operators over $m+1$ layers, like in convolution network architectures :

To reduce computations, wavelet convolutions are subsampled at intervals proportional to the last scale $2^{j_{q}}$, with an oversampling factor of $2$:

A final low-pass filtering and subsampling yields

With an FFT, the overall computational complexity is then $O(N\log N)$.

3 Scattering Distance and Deformation Stability

The scattering transform defines a distance between two images $f$ and $g$. This distance has important invariance and stability properties that are briefly reviewed. Let $|S_{J}f(x)|^{2}$ be the squared Euclidean norm of the vector $S_{J}f(x)$. The scattering distance of $f$ and $g$ is

For discrete images, the integral is replaced by a discrete sum. The scattering operator $S_{J}$ is contractive because it is a cascade of wavelet transforms $W_{J}$ and modulus operators, which are both contractive :

In particular $\|S_{J}f\|^{2}\leq\|f\|^{2}$. If the maximum order is $m=\infty$ then one can prove that if the wavelet transform is unitary then for appropriate complex wavelets $\|S_{J}f\|^{2}=\|f\|^{2}$. The energy of $f$ is thus spread across scattering coefficients of multiple orders, but this energy has a fast decay as the co-occurrence order $q$ increases. In the Caltech101 image database, 98% of the energy $\|S_{J}f\|^{2}$ is carried by scattering coefficients of order , $1$ and $2$. In applications, we shall thus limit the scattering order to $m=2$. The energy of all scattering coefficients of order $2$, $||f\star\psi_{j_{1},\gamma_{1}}|\star\psi_{j_{2},\gamma_{2}}|\star\phi_{J}$, is about 20% of the energy of all order $1$ coefficients $|f\star\psi_{j_{1},\gamma_{1}}|\star\phi_{J}$, which is not negligible. We shall see that order 2 coefficients have indeed an important impact on classification results.

The efficiency of a scattering representation comes from its invariance to local translations due to convolutions with $\phi_{J}$, and from its ability to linearize deformations. Let $D_{\tau}f(x)=f(x-\tau(x))$ be a deformation of $f$ with a regular displacement field $\tau(x)$. It is a pure translation only if $\nabla\tau=0$. We write $|\tau|_{\infty}=\sup_{x}|\tau(x)|$ the maximum translation amplitude, and $|\nabla\tau|_{\infty}=\sup_{x}|\nabla\tau(x)|$ the maximum deformation amplitude, where $|\nabla\tau(x)|$ is the matrix sup norm of $\nabla\tau(x)$. The sup-norm of the Hessian of $\tau$ is also written $|H\tau|_{\infty}$. It is shown in that the scattering metric satisfies

The first term $2^{-J}|\tau|_{\infty}$ is the translation error which is small if $2^{J}\gg|\tau|_{\infty}$. The other terms are dominated by the deformation amplitude $|\nabla\tau|_{\infty}$. If $2^{J}\geq|\tau|_{\infty}/|\nabla\tau|_{\infty}$ then two deformed signals have a scattering distance essentially proportional to the deformation amplitude $|\nabla\tau|_{\infty}$.

Classification by Affine Model Selection

A scattering representation $S_{J}f$ is invariant to small translations relatively to $2^{J}$. It linearizes deformations and provides co-occurence descriptors. A classifier is obtained by selecting an affine space model which best approximates $S_{J}f$.

Each signal class is represented by a random vector $F_{i}$ whose realizations are images of $N$ pixels in the class. Scattering vectors $S_{J}F_{i}(2^{J}n)$ define an image representation with a total of $N_{J}=2^{-2J}N\sum_{q=0}^{m}K^{q}\binom{J}{q}$ coefficients. Let $E\{S_{J}F_{i}(2^{J}n)\}$ be their expected values. Deformations of $F_{i}$ are mostly linearized by $S_{J}$ and thus produce a variability $S_{J}F_{i}-E\{S_{J}F_{i}\}$ which is well approximated in a linear space of low dimension $d$. This linear space is computed with a PCA by diagonalizing the covariance of $S_{J}F_{i}$. We denote by ${\mathbf{V}_{d,i}}$ the space generated by the $d$ covariance eigenvectors of largest variance. The dimension $d$ is adjusted so that $S_{J}F_{i}$ is closely approximated by its projection in the affine space

in comparison with the error produced by the affine spaces ${\bf A}_{d,i^{\prime}}$, $i^{\prime}\neq i$, corresponding to the other classes.

A signal $f$ will be associated to the class $\hat{\i}$ which yields the best affine space approximation:

where $\mathbf{V}_{d,i}^{\perp}$ is the orthogonal complement of $\mathbf{V}_{d,i}$. Minimizing the affine space approximation error is thus equivalent to minimize the distance between $S_{J}f$ and the class centroid $E\{S_{J}F_{i}\}$, without taking into account the first $d$ principal variability directions. A cross-validation procedure finds the dimension $d$ and the scale $2^{J}$ which yields the smallest classification error. This error is computed on a subset of the training images that is not used for the PCA calculations.

Affine space scattering models can be interpreted as generative models computed independently for each class. As opposed to discriminative classifiers such as an SVM, no interaction between classes is taken into account, besides the choice of the model dimensionality $d$.

Classification results are given for hand-written digits and textures that are deformed, rotated, scaled and have illumination variations. Scattering descriptors are computed with the complex Gabor wavelet (3) for $\xi=3\pi/4$, rotated along angles $k\pi/K$ with $0\leq k<K=6$. The lowpass filter is the Gaussian $\phi_{J}(x)=\lambda_{J}\exp(-(3x/2^{J+1})^{2}/2)$ with $\int\phi_{J}(x)dx=1$.

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. The state of the art is achieved with deep-learning convolutional networks and dictionary learning .

Table 1 compares the scattering PCA classifier at maximum orders $m=1$, $m=2$ and $m=3$. Cross validation finds an optimal scattering scale $2^{J}=2^{3}$. This value is compatible with observed deformations of digits whose amplitude is typically at most $8$ pixels. For $J=3$, there are $N/64$ second order scattering vectors $S_{J}f$ of dimension $127$ each.

Below $5\,10^{3}$ training samples, the scattering PCA classifier improves results of deep-learning convolutional networks. For $m=2$, second order scattering coefficients improve classification results obtained with $m=1$, but a third order $m=3$ scattering yields marginal improvements. An SVM classifier is also applied on scattering vectors for $m=2$, with a polynomial kernel whose degree was optimized. Minimum errors are obtained with a degree $4$. The SVM error is well above the PCA model selection error up to 60000 samples. For small training sets, it was indeed shown that generative models, which do not estimate cross terms between classes, can outperform discriminative classifiers such as SVM.

Table 2 gives the dimension $d$ of affine approximation spaces calculated by cross validation, for $m=2$. The normalized approximation error $\sigma_{d}^{2}$ is the expected approximation error $E\{\|S_{J}F_{i}-P_{\bf A_{i,d}}(S_{J}F_{i})\|^{2}\}$ in a class $i$ divided by the squared norm of $S_{J}F_{i}$, averaged over all $i$ and all $F_{i}$ in the test set. Table 2 shows that the cross-validation calculation of $d$ yields small approximation errors. Table 2 also gives the relative approximation error

produced by the closest affine model of a different class than that of $F_{i}$, averaged over all classes. As expected, when the training set increases, the dimension $d$ increases so $\sigma_{d}^{2}$ decreases and the relative approximation error $\lambda_{d}$ increases, which reduces the error rate.

Rotation invariance in the MNIST database is studied in the same setting as in . The authors have constructed a transformed database with 12000 training samples and 50000 test images, where samples are rotated versions of the digits using a uniform distribution in $[0,2\pi]$. The PCA incorporates rotation invariance by increasing the dimension $d$ of the affine space ${\bf A}_{i,d}$. It removes the main variability directions of $S_{J}f$ due to rotations. Error rates in Table 3 are smaller with a scattering PCA than with a convolution network . Better results are obtained with $m=2$ than with $m=1$ because second order coefficients maintain enough discriminability despite the removal of a larger number $d$ of principal directions.

The US-Postal Service dataset 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 4 gives results obtained with the PCA classifier and a polynomial kernel SVM classifier applied to scattering coefficients. The scattering scale was also set to $J=3$ by cross-validation.

2 Scattering Texture Classification

Scattering coefficients provide new texture descriptors, carrying co-occurrence information at different scales and orientations. A texture can be modeled as a realization of a stationary process $F(x)$. Scattering coefficients $S_{J}F(x)$ are obtained with successive convolutions and modulus operators which preserve stationarity. Averaging by $\phi_{J}$ does not modify expected values so $E\{S_{J}F(x)\}$ is a vector whose coefficients do not depend upon $x$ and $\phi_{J}$. The convolution with $\phi_{J}$ reduces the coefficient variability and for a large class of ergodic processes, the variance of $S_{J}F(x)$ decreases exponentially to zero as $J$ increases. As a result, $S_{J}F(x)$ is a good estimator of $E\{S_{J}F(x)\}$ when $J$ is sufficiently large. Figure 1 shows an example of such vector for a textured image with $m=3$.

Textures having same mean and same power spectrum have nearly the same scattering coefficients of order $q=0$ and $q=1$. However, different textures typically have co-occurence coefficients of order $q\geq 2$ which are different. Let $S_{q,J}F_{i}$ be the vector of scattering coefficients of order $q$ for a texture $F_{i}$. The distance of scattering vectors of order $q$ for two textures $F_{1}$ and $F_{2}$ is normalized by their variance $\sigma^{2}(S_{q,J}F_{i})$:

Table 5 gives $\rho_{q}(F_{1},F_{2})$ for two Brodatz textures in Figure 2, which have different power spectrum. Their expected scattering vectors $E\{S_{J}F_{q,i}\}$ have a relatively large distance $\rho_{q}(F_{1},F_{2})$ at all orders $q\geq 1$. The texture $\widetilde{F}_{1}$ in Figure 2 has same power spectrum as $F_{2}$. When $q=1$, equalizing the power spectrum reduces $\rho_{q}(\widetilde{F}_{1},F_{2})$ to (up to estimation errors) but $\rho_{q}(\widetilde{F}_{1},F_{2})$ remains well above zero for $q>1$. Textures having same power spectrum can thus be discriminated from scattering coefficients of order $q>1$.

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 it challenging for classification. Pose variations require global rotation invariance. Figure 3 illustrates the large intra class variability, and also shows that the variability across classes is not always important.

State of the art on this database achieves a 2.46% error rate, obtained in with an optimized Markov Random Field model. The scattering PCA classifier has a 0.09% error rate, which is a factor 25 improvement, as shown in Table 6. The database is randomly split into a training and a testing set, which either comprises 46 training images each as in , or contains 23 training images as in . Results are averaged over 10 different splits.

The cross-validation adjusts the scattering scale $2^{J}=2^{7}$ which is the maximum value. Indeed, these textures are fully stationary and increasing the scale reduces the variance of the scattering coefficients variability across realizations. Global invariance to rotation and illumination is provided by the PCA affine space models. They include the main variation directions of scattering vectors due to rotations or illumination variations.

The dimension of affine approximation space models is adjusted by cross validation to $d=6$ and $d=22$ respectively for $23$ and $46$ training samples. The resulting error rates are respectively $0.9\%$ and $0.09\%$. With an SVM using a polynomial kernel, the classification error for 46 training samples per class increases to $1.1\%$. The intra class normalized approximation error $\sigma^{2}_{d}$ is only $2.5\cdot 10^{-3}$ when using $46$ training samples, about half of the error produced in the case of $23$ training samples, in which $\sigma^{2}_{d}$ is $5.3\cdot 10^{-3}$. The estimated separation ratio is $\lambda_{d}=8$ and $\lambda_{d}=5$ respectively. Such low approximation errors are possible thanks to the fast variance decay of scattering coefficients as the scale increases and to the global invariance properties provided by the affine spaces.

Conclusion

A scattering transform provides a locally translation invariant representation, which linearizes small deformations, and provides co-occurrence coefficients which characterize textures. For handwritten digit recognition and texture discrimination with small training size sequences, a PCA model selection classifier yields state of the art results.

Besides translations, invariance can be extended to any compact Lie group $G$, by combining another scattering transform defined on $G$. The cascade of wavelet transforms in $\bf L^{2}(R^{2})$ is then replaced by a cascade of wavelet transforms in ${\bf L^{2}}(G)$ .

References