Optimally damped rank-reduction method

A shortcoming of the traditional rank-eduction and damped rank-reduction methods is that the output perfomance is quite sensitive to the input parameter $N$. The rank parameter varies greatly for different datasets, especially for field datasets. In order to alleviate the influence of the rank parameter, we take advantage of an adaptive singular-value weighting algorithm developed by Benaych-Georges and Nadakuditi (2012) and Nadakuditi (2013). An adaptive singular-value weighting matrix can be obtained by solving the following problem:

$$\hat{\boldsymbol{\omega}}=\arg\min_{\boldsymbol{\omega}} \paralle... ...i^S(\mathbf{v}_i^S)^H - \omega_i\sigma_i\mathbf{u}_i\mathbf{v}_i^H \parallel_F,$$ (11)

where $\mathbf{u}_i^S$ and $\mathbf{v}_i^S$ denote $i$th left and right singular vectors corresponding to the estimated signal component, $\sigma_i$ denotes the $i$th singular value in $\boldsymbol{\Sigma}$. $\boldsymbol{\omega}$ is the adaptive weight vector, which can be used to construct the adaptive weighting matrix $\hat{\mathbf{W}}=$diag$(\hat{\boldsymbol{\omega}})$. The solution to the optimization problem can be obtained as:

$$\hat{w}_i=\left(-\frac{2}{\sigma_i}\frac{\mathcal{D}(\sigma_i;\boldsymbol{\Sigma})}{\mathcal{D}'(\sigma_i;\boldsymbol{\Sigma})}\right).$$ (12)

where $\mathcal{D}$ represents a transform expressed as:

$$\mathcal{D}(\sigma;\boldsymbol{\Sigma})= \frac{1}{N^2} \left(\sigma(\sigma^2\mathbf{I}-\boldsymbol{\Sigma}\boldsymbol{\Sigma}^H)^{-1}\right) \left(\sigma(\sigma^2\mathbf{I}-\boldsymbol{\Sigma}^H\boldsymbol{\Sigma})^{-1}\right),$$ (13)

and $\mathcal{D}'$ denotes its derivative. The expression of $\mathcal{D}'$ can be expressed as:

$$\begin{split} D'(\sigma;\boldsymbol{\Sigma})&= 2\left[\frac{1... ...athbf{I}-\boldsymbol{\Sigma}^2)^{-2}\right)\right]. \end{split}$$ (14)

The symbol Tr$(\cdot)$ denotes the trace of the input. The trace of a square matrix $\mathbf{X}$ is defined to be the sum of elements on the main diagonal of $\mathbf{X}$.

$$(\mathbf{X}) = \sum_{i=1}^{N}X_{i,i},$$ (15)

where $X_{i,i}$ denotes the main diagonal elements of $\mathbf{X}$.

The adaptive weighting operator can be applied to make the traditional rank-reduction method adaptive:

$$\hat{\mathbf{S}}=\mathbf{U}_N\hat{\mathbf{W}}\boldsymbol{\Sigma}_N\mathbf{V}_N^H.$$ (16)

To incorporate the adaptive weighting operator into the damped rank-reduction framework, we can introduce intermediate variables as

$$\mathbf{U}_N^{Q} = \mathbf{U}_N,$$ (17)

$$\boldsymbol{\Sigma}_N^{Q} = \hat{\mathbf{W}}\boldsymbol{\Sigma}_N,$$ (18)

$$\mathbf{V}_N^{Q} = \mathbf{V}_N,$$ (19)

then equation 16 turns into

$$\hat{\mathbf{S}}=\mathbf{U}_N^{Q}\boldsymbol{\Sigma}_N^{Q}\mathbf{V}_N^{Q}.$$ (20)

It can be derived that the damping formula (equation 8) also holds for equation 20 but with the damping threshold matrix expressed as

$$\boldsymbol{\Gamma} = \hat{\sigma}^K\left(\boldsymbol{\Sigma_N^Q}\right)^{-K},$$ (21)

where the subscript $N$ denotes a sufficiently large rank parameter, as required by the derivations detailed in Chen et al. (2019b). The resulted final form of the optimally damped rank-reduction method then can be expressed as

$$\hat{\mathbf{S}} = \mathbf{U}_N\mathbf{P}\hat{\mathbf{W}}\boldsymbol{\Sigma}_N\mathbf{V}_N^H.$$ (22)

The algorithm workflow for the optimally damped rank-reduction method (ODRR) is outlined in Algorithm 3. Due to the existence of a weighting matrix in the proposed method, it is more convenient to choose the input rank parameter $N$ empirically in practice. As mentioned previously, the rank parameter $N$ is in practice set to a relatively large value to forestall damaging the signal. But in the proposed algorithm, even if a large $N$ is used, the algorithm can adaptively shrink the singular-values. Thus, its performance is not sensitive to the input rank parameter, in contrast to other related approaches (Oropeza and Sacchi, 2011). Because of the insensitivity, one can use a sufficiently large $N$ in processing complicated datasets without leaving strong residual noise in the result. The convenience in tuning parameters makes the rank-reduction related methods more computationally feasible in large-scale data processing, e.g., the 5D reconstruction problem (Chen et al., 2019b), since one no longer needs to tune the parameters many times while one trial is already computationally demanding. A glossary describing the main mathematical notations are presented in Table 1. The three methods are all variations of how the singular-values are thresholded. Simply speaking, the DRR method improves the RR method by introducing a damping operation and the ODRR method improves DRR method by further introducing a weighting operation.

Table 1: Glossary of main mathematical notations used in this paper.

Symbols	Meanings
$\mathbf{D}(t/w,x,y)$	3D noisy seismic data
$\mathbf{D}(w_0)$ or $\mathbf{D}$	abbreviated notation of 2D frequency slice
$\mathbf{R}_i$	$i$th Hankel matrix
$\mathbf{M}$	block Hankel matrix
$\mathbf{S}$	signal component
$\hat{\mathbf{S}}$	estimated signal component
$\mathbf{U}$	left singular vector matrix
$\boldsymbol{\Sigma}$	singular-value matrix
$\mathbf{V}$	right singular vector matrix
$\hat{\mathbf{W}}$	weighting matrix
$\mathbf{P}$	damping matrix
$\boldsymbol{\Gamma}$	damping threshold matrix
$\mathbf{Q}$	temporary variable
$\mathcal{D}$	D-transform
$\mathcal{H}$	Hankelization operator
$\mathcal{T}$	thresholding operator
$\mathcal{A}$	averaging operator
$N$	rank
$K$	damping factor

Table 2: Comparison of SNRs in dB for different rank-reduction based methods. RR denotes the rank-reduction method. DRR denotes the damped rank-reduction method. ODRR denotes the optimally damped rank-reduction method.

Tests	Noisy (dB)	RR (dB)	DRR (dB)	ODRR (dB)
Linear synthetic (N=3)	-8.39	6.57	10.29	11.27
Linear synthetic (N=6)	-8.39	3.45	8.35	10.86
Hyperbolic synthetic (N=10)	-2.17	8.27	9.58	9.65
Hyperbolic synthetic (N=20)	-2.17	7.04	10.08	11.00

Table 3: Comparison of computational costs in seconds for different rank-reduction based methods.

Tests	FK (s)	RR (s)	DRR (s)	ODRR (s)
Linear synthetic (N=3)	0.13	2.27	2.45	2.43
Linear synthetic (N=6)	0.17	3.02	3.04	3.17
Hyperbolic synthetic (N=10)	0.43	1.69	1.71	1.83
Hyperbolic synthetic (N=20)	0.49	2.75	2.69	2.83