Seislet-based morphological component analysis using scale-dependent exponential shrinkage

Sparsifying MCA with exponential shrinkage shaping

The IST algorithm used by MCA requires soft thresholding function to filter out the unwanted small values. Besides soft thresholding (Donoho, 1995), many other shrinkage functions can also be applied to obtain possibly better sparseness. One particular choice is hard thresholding:

$$t_{\lambda}(u)=\mathrm{Hard}_{\lambda}(u)=u.*(\vert u\vert>\lambda?1:0).$$ (19)

where $(\cdot)?A:B$ frequency used hereafter is an if-else judgment in C-code style: The expression equals $A$ if the statement $(\cdot)$ is true, and $B$ otherwise.

Another choice is Stein thresholding (Peyre, 2010; Mallat, 2009):

$$t_{\lambda}(u)=\mathrm{Stein}_{\lambda}(u)=u.*\max\left(1-(\frac{\lambda}{\vert u\vert})^{2},0\right).$$ (20)

Stein thresholding does not suffer from the bias of soft thresholding, that is,

$$\vert\mathrm{Stein}_{\lambda}(u)-u\vert \rightarrow 0, \vert\mathrm{Soft}_{\lambda}(u)-u\vert\rightarrow \lambda, \mathrm{if}\; u\rightarrow\infty.$$ (21)

Recent advances in nonconvex optimization (Chartrand, 2012; Voronin and Chartrand, 2013; Chartrand and Wohlberg, 2013) show that the shrinkage operator in IST algorithm (Eq. (2)) can be generalized to a $p$ -quasinorm ($0<p\leq1$ ) thresholding operator $T_{\lambda}$ , in which

$$t_{\lambda}(u)=\mathrm{pThresh}_{\lambda,p}(u)=u.*\max\left(1-(\frac{\lambda}{\vert u\vert})^{2-p},0\right).$$ (22)

A special case is that of $p=1$ , which corresponds to the soft thresholding operator exactly.

Most of these shrinkage functions interpolate between the hard and soft thresholders. It is tempting for us to design a more general shrinkage function to sparsify the transform domain coefficients in shaping regularized MCA. One possibility is multiplying an exponential factor on the elements of original data:

$$t_{\lambda}(u)=u.*\mathrm{exp}(-(\frac{\lambda}{\vert u\vert})^{2-p}).$$ (23)

Based on Taylor series, this operator in Eq. (23) enjoys some useful properties:

• It is valuable to point out that the exponential shrinkage can be considered as a smooth $\ell_0$ constraint (Gholami and Hosseini, 2011; Mohimani et al., 2009). For $\vert u\vert>>\lambda$ , it is a good approximation of the p-thresholding operator in Eq. (22), and does not suffer the bias when $p\neq 1$ . It reduces to Stein thresholding operator for $p=0$ and soft thresholding for $p=1$ .

  $$\begin{split}t_{\lambda,0}(u)=u.*\mathrm{exp}(-(\frac{\lambda... ...rt}))\approx u.*(1-(\frac{\lambda}{\vert u\vert})). \end{split}$$ (24)

  
  

• It is free of non-differentiable singularity at the thresholding point $\lambda$ . The transition between the small values and the large values is smoothly stretched. Due to the exponential factor less than 1 ( $\mathrm{exp}(-(\frac{\lambda}{\vert u\vert})^{2-p})<1$ ), this operator will slightly decrease the data amplitude, even for $\vert u\vert<\lambda$ .
• Besides the threshold $\lambda$ , we have another independent parameter $p$ which can be flexibly chosen to achieve better performance.

In the language of shaping regularization, shrinkage-based shaping operator $S$ is equivalent to multiplying the coefficient vector $x$ by a diagonal weighting matrix $W$ to in the sense that

$$S(x)=Wx,$$ (25)

where

$$\mathrm{diag}(W_{ii})= \begin{cases}1-\frac{\lambda}{\vert x_i\ve... ...xp}(-(\frac{\lambda}{\vert x_i\vert})^{2-p}),& \mathrm{Exponential} \end{cases}$$ (26)

For the convenience of comparison, we plot these thresholding operators in Figure 1.

Note that we are using seislet transform which has different scales for signal representation. Usually, large scales of seislet coefficients corresponds to unpredictable noise, while most of the important information gets transformed into smaller scales. We design a scale-dependent diagonal weighting operator:

$$\mathrm{diag}(W_{ii})= s(x_i)<s_0?1:0,$$ (27)

where $s_0$ is user-defined scale, while $s(x_i)$ is the scale that the coefficient $x_i$ correspond to. Putting all things together, in the MCA shaping regularization we are using a scale-dependent exponential shrinkage operator which is a composite operator cascaded with a scale-muting operator $W_{s}$ and an exponential weighting operator $W_{exp}$ :

$$S(x)=W_{exp} W_s x.$$ (28)

The use of scale-dependent exponential shrinkage offers easy control on the separation of the signal components we would like to capture. It is interesting to mention that under the Fourier basis, the scale-muting operator $W_s$ becomes a frequency mask (it behaves like a selective hard thresholding), which can be employed to remove the groundroll in application to seismic interpolation (Gholami, 2014).

By incorporating PWD dip estimation and scale-dependent exponential shrinkage shaping, we summarize the proposed seislet-based MCA algorithm as Algorithm MCA. Seislet transforms associated with different dips form a combined seislet frame (Fomel and Liu, 2010). The threshold in each iteration can be determined with a predefined percentile according to Hoare's algorithm. Shrinkage operator plays the role of crosstalk removal in MCA algorithm, as explained in more detail in Appendix A.

thresh
Figure 1. A schematic plot of the shrinkage operators,$\lambda =1$

Seislet-based morphological component analysis using scale-dependent exponential shrinkage