Wavelet Scattering properties & implementation?

Question

What are the properties of the scattering transform, besides shift invariance and warp stability? How is it implemented in practice, and how can one visualize its computational graph?

Answer 1

Scattering overview provided in this answer.

Computational structure

_{Fig 4, Deep Scattering Spectrum}

In steps:

1. (First order begins) $x$ convolves with $\psi1_i$ --> $W1_i$
2. Modulus, $W1_i \rightarrow U1_i$
3. Lowpass, $U1_i \rightarrow S1_i$
4. Repeat for all $i$
5. (Second order begins) $U1_0$ convolves with $\psi2_j$ --> $W2_{0_j}$
6. Modulus, $W2_{0_j} \rightarrow U2_{0_j}$
7. Lowpass, $U2_{0_j} \rightarrow S2_{0_j}$
8. (Second order continues) $U1_1$ convolves with $\psi2_j$ --> ... --> $S2_{0_j}$
9. (Second order continues) Repeat for all $i, j$
10. (Third order begins) $U2_{i_j}$ convolves with $\psi3_k$ --> ... --> $S3_{i_{j_k}}$
11. (Third order continues) Repeat for all $i, j, k$
12. ...

That is, at each order $m$, we convolve every coefficient from previous order with every wavelet of current order; the number of coefficients thus grows exponentially. Zeroth-order discussed in "Energy analysis".

The idea's simple: at first order, we extract frequencies: oscillation rates of the signal. At second order, we extract rates of those rates; to do so, we treat every rate from first order as its own signal, and repeat first order upon each. Except, it's not upon CWT directly, but its modulus, i.e. amplitude modulations - so second order is rates of AMs.

Key differences with CNNs: 1) kernels are fixed, not learned; 2) every layer's (order's) coefficients comprise the "output" (black dots in diagram).

Pseudocode

Unoptimized implementation is fairly straightforward:

U1, S1, S2 = [], [], []

# zeroth order
S0 = conv(x, phi_t)  # (1, time)

# first order
for p1t in psi1_t:
    u1 = abs(conv(x, p1t))
    S1.append(conv(u1, phi_t))
    U1.append(u1)
S1 = array(S1)  # (freq, time)
U1 = array(U1)

# second order
for p2t in psi2_t:
    u2 = abs(conv(x, U1))
    S2.append(conv(u2, phi_t))
S2 = array(S2)  # (AM freq, freq, time)

Increased frequency resolution

For the same mean frequency resolution, CWT is worse than STFT at high frequencies. Scattering increases the effective resolution by extracting, via higher orders, the differences of frequencies: if we cannot distinguish $\xi_2$ from $\xi_1$, we may still find $|\xi_2 - \xi_1|$ (Ref 1, VI.B).

Borrowing an example from Ref 1, Fig 6:

• First plot is scalogram, second is first-order scattering, third is second-order.
• Left signal at $\xi_1$ is a $600 \text{ Hz}$ and $675 \text{ Hz}$ chord (played simultaneously), right same frequencies but appregio (played in quick succession)
• Scalogram shows the chords blended together, and $S_1$ makes chords and appregios look the same per time averaging
• $S_2$ resolves the blended frequencies

This resolving is owed to the fact that, closely spaced tones are equivalently a single amplitude-modulated tone at their mean frequency:

$$ \cos(A) + \cos(B) = 2 \cos((A - B)/2) \cos((A + B)/2) \tag{1} $$

where $(A + B)/2$ is the "carrier" and $(A - B)/2$ the "modulator". Thus if $A$ and $B$ (or $\xi_2$ and $\xi_1$) are close, then we have an A.M. whose envelope varies with frequency $|\xi_2 - \xi_1|/2$ (the half dropped in plot), and the $S_2$ will extract it just like $S_1$ would extract $|\xi_2 - \xi_1|$ from $x(t) = \cos(|\xi_2 - \xi_1|t)$. (This is an approximate picture which is sometimes exact; for full description see Ref 1, VI B & C).

Subsampling

Is lossess "within tolerance", since wavelets aren't perfectly compact. Though, they can be: Morlets decay per $e^{-\omega^2}$, which hit sub-machine epsilon quickly. Thus convolving sampled can equal sampling a (continuous-time) convolution result, within float precision.

In practice a reasonable threshold is picked; Kymatio uses criterion_amplitude=1e-3, meaning, the wavelet's "support" is defined to end at the first sample which drops below 1/1000 of peak amplitude, or $10^{-6}$ of peak energy. We can subsample by any factor, but work only with powers of 2 for convenience. Taking the previous filterbank:

• blue = can't subsample without aliasing; orange = can subsample by $2^1$; green = by $2^2$
• vertical black lines = octave separators
• note lowpass is contained left of fourth line, meaning $2^3 = 8$. This is $T$. In implementation, $J$ and $T$ are usually global, so all orders share the lowpass. Lowpassing by $\phi_T$ thus allows subsampling all coefficients by $T$, and concatenating to form an array feedable to e.g. an NN

Boundary effects

Using zero padding and boxcar to show filter "support" for clarity:

Minimal ideas with avoiding boundary effects are:

1. Left shall not draw from right: don't want beginning described by what happened 10 mins later
2. Wavelets fully decay: large scale wavelets' support can easily exceed input's length
3. Convolutions are energy-preserving (optional): zero padding wrecks this for large scales, as wavelets multiply with zeros

One may find the $T$ parameter surprising visually:

The support of the lowpass, defined by criterion_amplitude=1e-3, is x16 the width (i.e. scale), $T$. If we zoom on the $T$ portion, it looks nearly flat. Same for wavelets, largest scale being $2^J$. (Note, this is temporal support, relevant for boundary effects, whereas for subsampling we use frequential support)

Padding

Besides avoiding boundary effects, the idea's to emulate no padding - i.e. having an infinite signal (rather, longer than largest-scale wavelet's). It comes down to picking the "best guess" for the signal's continuation, which is best done with knowledge of the source process.

There isn't a one-fits-all, but my recommendation for the general case is reflect padding, which extends the signal from its local features, and is energy-conserving.

Energy analysis

Flow of energy is flow of information. Write the filterbank as${}^{1}$

$$ A(\omega) = |\hat\phi (\omega)|^2 + \sum_{\lambda \in \Lambda}|\hat\psi_\lambda(\omega)|^2,\ \forall \omega \in \mathbb{R} \tag{2} $$

where $\Lambda$ is the set of all center frequencies. This is Littlewood-Paley summation, the sum of all filters' energies, which satisfies

$$ 1 - \alpha \leq A(\omega) \leq 1,\ \alpha < 1. \tag{3} $$

Complete wavelet transform writes

$$ Wx = \left(x \star \phi, \vphantom{\frac{.}{.}} x \star \psi_\lambda\right)_{t\in R, \lambda \in \Lambda} \tag{4} $$

that includes the lowpass to capture DC, which in a finite implementation requires infinite scale. Applying Parseval-Plancherel's theorem:

$$ (1 - \alpha) ||x||^2 \leq ||Wx||^2 \leq ||x||^2 \tag{5} $$

where $||x||^2 = \int|x(t)|^2dt$ is signal energy, and $||Wx||^2 = \int |x \star \phi|^2 dt + \sum_{\lambda \in \Lambda}\int|x \star \psi_\lambda|^2 dt$ is sum of coefficient energies. If $\alpha=0$, then the filterbank is a tight frame: input energy is conserved exactly.

• The upper bound of $\leq1$ guarantees the filterbank is contractive: output energy doesn't exceed input's (else cascaded filtering could blow up energy). It also conforms to scattering being a decomposition rather than generative.
• The lower bound of $>0$ ensures $Wx$ has an inverse (else there exists some input for which the transform yields nothing, and we can't invert nothing). $\gg 0$ additionally ensures the inverse is numerically stable.

Taking an example filterbank and its LP-sum,

Exceeding 1 is an addressable implementation artifact, and so are sum oscillations. Assume for simplicity we have $\alpha = 0$; then the lowpass-wavelets energy split is

Usually the black portion (lowpass) is much lesser, but for clarity suppose it's even greater. Then, across scattering orders, the energy flow looks like:

where left shows fixed proportions for all orders, and right shows "realistic" proportions. Verbally, each order: 1) takes wavelet transform of previous order, takes its modulus, and passes to next order; 2) outputs lowpassing of previous order. $S_2$ is thus a second-order coefficient, produced by second-order's lowpassing of first-order's bandpassing of zeroth-order coefficients (i.e. the signal).

• Right is realistic because the complex modulus shifts all coefficient frequencies to be globally lower (Ref 7). In practice the DC is most dominant, followed by neighboring frequencies. Empirical results in Ref 1, Table 1.
• $T$ controls the proportions: greater $T$ reduces energy for each order's output, $S$, and increases for $U$. Intuitively, it takes more orders to recover information, as more information is lost to averaging. Further, by design, i.e. Eq $(2)$, raising $T$ shrinks lowpass's frequency support and hence its energy slice.

1: here $\hat\psi_\lambda$ is shorthand for analytic and anti-analytic wavelets, depending on whether $\lambda$ is positive or negative. In implementation, for real inputs, we do only analytic and normalize the sum to 2 for just analytic, so the sum reads $|\hat\phi(\omega)|^2 + \frac{1}{2}\sum_{\lambda \in \Lambda} \left(|\hat\psi_\lambda(\omega)|^2 + |\hat\psi_\lambda(-\omega)|^2 \right)$.

Analysis vs Synthesis, 2.0

Flow of energy is flow of information: scattering has infinite layers, yet most implementations stop at two, citing that higher orders have "negligible energy". So what if they do?

Suppose we take the complete scattering transform, and $||S_0|| + ||S_1|| + ||S_2|| = .999||x||$. Orders $3$ to $\infty$ hence contain 0.1% of input's total energy. Next suppose we attempt to recover $x$ by inverting all scattering coefficients: perfection (within global phase). Now, suppose we discard all orders $3$ to $\infty$ and invert only up to $S_2$: the result will differ in energy by at most 0.1%.

That is, $S_3$ to $S_\infty$ describe only a small fraction of the input; it's as if we applied a transform on only this tiny "leftover". Could it possibly contain important information? Sure - but unlikely. The practical ratio for first two orders is also much less, (.7, .8, .9, depending on signal and $T$; Ref 1). Thus, analysis and synthesis information are fundamentally tied.

Note: "analysis vs synthesis" is also a distinction in frame theory, relevant to scattering. The word use here is my own and differs from it (except "synthesis" in both contexts means recovery).

Invertibility

Myth: scalogram isn't invertible because modulus kills phase. This is true of the Fourier transform, which is orthogonal. The CWT is a highly redundant transform, and phase information persists even past modulus.

Invertibility and its criterion were proven by S. Mallat and I. Waldspurger. We know a scalogram cannot be perfectly invertible due to phase-shift invariance; instead, it is invertible within global phase shift, $e^{j\omega_0}$.

This is a strong general invertibility. A mono-component signal, or signal writable in form $A(t) \cos(\omega t)$, is invertible within a simple phase shift:

1. $A(t)$ is recovered perfectly directly from scalogram
2. $\phi(t)$ is recovered within global phase: $\tilde{\phi}(t) = \phi(t) + C$
3. We recover $A(t)e^{j\tilde{\phi}(t)}$, whose real part is $\boxed{A(t)\cos(\phi(t) + C)}$, where $A(t)^{j\phi (t)}$ is the analytic part of $x(t)$

Thus, instantaneous amplitude and frequency are recovered perfectly. In most cases, we don't care for such phase shift. If the signal is multi-component and its components can be fully separated, which holds sufficiently for a wide variety of signals, this invertibility holds. Else, the most general signal may degrade severely by time-domain appearance - but time-frequency ridges, and amplitude, will still recover perfectly. Complementary methods (e.g. synchrosqueezing) can improve component separation and relative phase recovery.

Scattering is invertible within a global phase shift and a global time shift $T$.

Reconstruction: gradient-based

Inverting scalogram of exponential chirp:

Can do same with scattering, and we'll observe worse reconstruction with greater $T$ unless we include higher orders; this is illustrated directly on audio in S. Mallat's lecture (17:32).

Full post

Time-warp equivariance

We find insight comparing impulse response${}^1$ of CWT and scattering, latter for all $T$; peak amplitude of each row has been normalized to be same for visual clarity (hence CWT', S'):

^{(1: unlike usual, scattering's IR can't be used as surrogate in convolution (i.e. overlap-adding weighted IRs) due to modulus nonlinearity)}

Elaborated in this post

Filterbank implementation

In addition to all discussed, worth noting for the kymatio implementation:

1. $J$ controls number of octaves, and $Q$ the number of wavelets per octave (which indirectly sets the actual quality factor). The total number of wavelets is thus $J \cdot Q$.
2. The lowest frequencies (typically last octave) are tiled by STFT rather than CQT - that is, wavelets with fixed bandwidth and linear center frequency distribution. This avoids discretization artifacts and enables complete tiling with a finite number of filters. Coupled with other quality checks, the total number of wavelets may differ significantly from $J \cdot Q$, especially for low $J$ (significant only for undershooting - excess is minimal and rare)
3. Only the analytic (left-half) part of the filterbank is sampled, per assuming input is real-valued; if input is complex, the first order will include the anti-analytic part, but not higher orders (modulus is always real).
4. Morlet is pseudo-analytic; it leaks into negatives near Nyquist and dc. However, the negatives can be forced to zero - at expense of nicer temporal decay, which is often worth it.

Example: Amplitude Modulated Cosine

Filterbanks

Signal

CWT, first order

CWT, second order

First order correctly identifies the carrier at $66.8 \approx 64\ \text{Hz}$ (closer with greater $Q_1$ that'd tighter tile the frequency), and second order the modulator at $3.2 \approx 3\ \text{Hz}$ (closer with greater $Q_2$).

Example: White Gaussian Noise

Filterbanks

Same

Signal

CWT, first order

CWT, second order

Energy analysis (cont'd)

Note, highest xi2 are mostly "empty", with only top rows for each xi2 (high xi1) having notable energy. This is consistent with the global energy shifting toward lower described in Ref 7. Intuitively, a $5\ \text{Hz}$ carrier cannot have a $50\ \text{Hz}$ amplitude modulator (else the modulator is actually the carrier), nor can a $50\ \text{Hz}$ carrier have a $45\ \text{Hz}$ modulator as that gets resolved as separate frequencies (Eq $(5)$). Per theorems, they should not exist at all (but do because of discretization and finiteness).

Thus frequencies near Nyquist (xi ~ .5) will lack AMs that are near Nyquist, and those not near it. In fact the highest xi2 shown above is .2, which is far from Nyquist; the actual highest was omitted (since 10 is awkward to subplot).

Implementations account this via a conditional: if xi2 < xi1: compute. This is akin to imposing $f_\text{modulator} < f_\text{carrier}$. Nonetheless, non-AM waveforms can have high xi2 contents if first-order wavelets are highly time-localized (e.g. square wave), but in most configurations xi2 < xi1 is a safe assumption.

Properties summary

Some are explained in Ref 1, others in these answers. Division into "reducing/decreasing/linearizing" isn't strict. Most properties are proven for both deterministic and stochastic processes.

Variability-reducing

1. Time-shift invariance: $S(x(t - c)) \approx S(x(t)),\ c \ll T$
2. Contraction: $||Sx - Sy|| \leq ||x - y||$
3. Energy preservation: $||Sx|| = ||x||$, for tight frame filterbank. Energy decay (hence convergence) is exponential across orders (Ref 7)
4. Phase-shift invariance: $S(x(t)) = S(y(t)),\ y(t) = \mathcal{F}^{-1}\left\{e^{-j\omega_0} \hat{x}(\omega)\right\}$ (due to modulus)
5. Noise-robustness (additive): follows from 2
6. Invariance to filtering with regularities: $\tilde S_2 (x \star h)(t, \lambda_1, \lambda_2) \approx \tilde S_2 (x)(t, \lambda_1, \lambda_2)$, if $\hat h(\omega)$ is approximately constant on the support of $\hat\psi_{\lambda_1}$. Ref 1, VI.A.
   • $\tilde S$ is normalized scattering: $\tilde S_2 = S_2 / (S_1 + \epsilon)$
   • "regularities" refers to the "approx constant" criterion - which is, precisely, $\lambda_1/Q_1 \ll \left(\int |t||h(t)|dt\right)^{-1}$

Variability-increasing

1. Amplitude modulations: $A(t)$ in $A(t)\cos(\omega(t))$
2. Amplitude modulation rates: $A'(t)$
3. Frequency modulations: $\omega(t)$
4. Frequency modulation differentials: $|\omega_1(t) - \omega_0(t)|$

• "increasing" as in supplying discriminative information
• 2 and 4 propagate through every order: 3rd order captures rates of amplitude modulation rates, and differences of differences of frequency bands; 4th captures rates of rates of rates, etc.

Variability-preserving

1. Translation equivariance: $x(t) \rightarrow x(t - c) \Leftrightarrow S(t) \rightarrow S(t - c)$. Invariance holds for $c \ll T$ in L2 and $c < T$ in subsampling sense, but a global shift still shifts the entire representation and non-$(c \ll T)$ structures.
2. Scale equivariance: $\left(\mathcal{W}_{(a,b)}x(t) \star \psi\right)(t) = \mathcal{W}_{(a,b)}\left(x(t) \star \mathcal{W}_{(1/a,0)}\psi\right)(t)$, for structures $\gg T$ in support (but holds globally for CWT) -- scaling input (in amplitude and time) corresponds to scaling its time-frequency representation. Explained here (coming soon)
   • $\mathcal{W}_{(a,b)}x(t) = ax(at + b)$
3. Time-warp equivariance (multiplicative): $(\hat{x}(\omega))' \rightarrow (\hat{x}(\omega/c))' \Leftrightarrow S'(t, \lambda) \rightarrow S'(t, |\lambda| - \log(c)),$ for $a_\lambda < \text{log}(T)$. Warping a signal has the ~same effect on its time-frequency representation for any signal's location on time-frequency plane.
   • Equivalently, $S(x(\tau(t) t) - S(x(t)) \approx S(y(\tau(t) t) - S(y(t)),\ y(t) = x((1 - \epsilon)t)$ ($y$ is constant multiplicative warping of $x$; $\epsilon$ is such that $(1 - \epsilon) > 0$ - doesn't need to be small, can be negative)
   • $x' = x(\tau(t)t)$, and $a_\lambda=$ scale of log-frequency $\lambda$
   • $a_\lambda < \log(T)$ is always satisfied in standard scattering per $\log(T) = J$
   • $(\hat{x}(\omega / c))'$ means first warping by $c$ then by $\tau(t)$; former is a simple frequency scaling
   • This is a result I derived, and haven't seen a theorem for - my argument here (coming soon).
4. (Unverified) Time-warp equivariance (additive): 3 but $(\tau(t) t) \rightarrow (t - \tau(t))$, with $\max|\tau'| \ll 1/Q$ and $\max|\tau| \ll T$. Confirmed with some numeric experimentation but I've not verified further.
5. Ridges, modulations, other: separable representations in time-frequency for some classes of signals. Further linearity is obtainable via additional nonlinear processing steps (e.g. learned convolutions, JTFS).

Variability-linearizing

1. Time-warp stability: $||S(x(t - \tau(t))) - S(x(t))|| \leq CK(\tau(t)) ||x||,\ |\tau'(t)| < 1$
   • $C=$ constant determined by filterbank, $\propto Q$
   • $Q = $ (center frequency) / bandwidth -- Constant Q Transform
   • $K=$ a function of $\tau$ involving derivatives and scale of invariance (Ref 3 Eq 50)
   • Linearizing holds for $\max|\tau| \ll T$ and $\max|\tau'| \ll 1/Q$
   • The theorem assumption may be $|\tau'| < 1/Q$ instead, or other dependence on $Q$, currently unclear; my experiments verify $<1$. Will update

Other

1. Sparsity: most energy is concentrated over few coefficients for a wide variety of signals, even if noised
2. Spatial coherence: coefficients can be convolved over in time, in log-frequency, and jointly, enabling fixed and learnable higher-order modulation and invariant extractions
3. Unmixing of offsets & modulations: unlike STFT, CWT enforces zero-mean kernels, keeping modulations separate from (DC) offsets, with important features implications. Equivalently, unmixing rates and shifts, or weights and biases.
4. Configurability: scales of invariance, linearity, and stability, and redundancy and time-frequency tradeoff of representation, can be adapted to application. Invariance can also be learned.
5. Differentiability: all involved operators are differentiable, hence scattering can be inserted as an intermediate block in an NN
6. Dimensional generality: scattering preserves its invariance and discriminability properties for higher dimensional tasks (images, volumes, etc), and even gains new.

"Linearizing"

Formally, means that representations of deformed signals form a ~linear subspace of undeformed. This can be measured with a discretized Laplacian:

$$ \nabla^2_\tau Sx(t, p) = Sx(t, p) - \frac{1}{2} \left(Sx_{\tau}(t, p) + Sx_{\tau^{-1}}(t, p) \vphantom{\frac{.}{.}}\right) $$

where $\nabla^2_\tau$ differentiates along the deformation, and $x_\tau(u) = x(u - \tau(u))$.

• Interprets as, equal and opposite "nudges" of $\tau$ average into un-nudged in representation
• Implies all such deformations can linearly combine into undeformed, i.e. are scalar-scaled versions of each other
• Distinct from equivariance and invariance, but possible to confuse with latter if not careful as "nudges" yield small distance by definition

"Invariance" vs "Stability"

"Invariance" means coefficient distance remains small, and changes little, for any deformation size under the scale of invariance. In theoretic work it often means this distance is and stays zero.

"Stability" means coefficient distance is bounded by only by size of deformation, not characteristics of the signal (besides its norm). That is, we can't find $x$ whose coeff distance with $x'$ exceeds said bound; this isn't true of STFT, where we can always drive distance to maximum ($||x|| + ||x'||$ per contraction) by frequency-shifting $x$, nor of CWT per time-shifting.

• Both formulate via Euclidean distance, $||\Phi(x) - \Phi(x')||$, but "stability" includes the size of deformation in its bound, whereas scale of invariance doesn't (it's replaced by an upper bound on deformation size).
• Stability doesn't reduce variability, so discriminative power isn't lost. Instead, it bounds variability, and for small deformations, linearizes. Discriminability is increased.
• Stability is related to warp equivariance: coefficients change ~same for same deformation, regardless of signal's characteristics.
• "Size" of deformation includes its norm and shape (derivatives).
• "Stability", formally, means Lipschitz continuity.

T Note: statements on scattering being "invariant over the scale of invariance", are not correct, unless interpreted in a specific way not stated here or in any reference I know of. _{"Scale of invariance" is not a thing in the sense described, nor can it be: a window is smooth and continuous without a strict cutoff, and said invariance would in fact contradict the motivation of the stability theorem (a shift is a warp). The statements follow my discussions with Kymatio, but in fact are absent from the founding paper, Graph Invariant Scattering. I may elaborate later.}

Conclusion

Scattering is dope