A Domain Decomposition Fourier Continuation Method for Enhanced $L_1$ Regularization Using Sparsity of Edges in Reconstructing Fourier Data

Abstract

$L_1$ regularization is widely used in various applications for sparsifying transform. In Wasserman et al. (J Sci Comput 65(2):533–552, 2015) the reconstruction of Fourier data with $L_1$ minimization using sparsity of edges was proposed—the sparse PA method. With the sparse PA method, the given Fourier data are reconstructed on a uniform grid through the convex optimization based on the $L_1$ regularization of the jump function. In this paper, based on the method proposed by Wasserman et al. (J Sci Comput 65(2):533–552, 2015) we propose to use the domain decomposition method to further enhance the quality of the sparse PA method. The main motivation of this paper is to minimize the global effect of strong edges in $L_1$ regularization that the reconstructed function near weak edges does not benefit from the sparse PA method. For this, we split the given domain into several subdomains and apply $L_1$ regularization in each subdomain separately. The split function is not necessarily periodic, so we adopt the Fourier continuation method in each subdomain to find the Fourier coefficients defined in the subdomain that are consistent to the given global Fourier data. The numerical results show that the proposed domain decomposition method yields sharp reconstructions near both strong and weak edges. The proposed method is suitable when the reconstruction is required only locally.

Appendix

The key element of the FC method is to find a smooth matching function. The matching function is not unique but depends on various parameters used. The matching function $f_{\text {match}}(z)$ is approximated with the Gram polynomials for conditioning purpose defined as below. For more details, see [15].

Definition 2

Gram Polynomials: The Gram polynomials $P_n(\xi )$, where $-1\le \xi \le 1$, satisfy the following three-term recurrence relation,

$$\begin{aligned} P_n(\xi )= & {} \alpha _{n-1} \xi P_{n-1}(\xi ) -\frac{\alpha _{n-1}}{\alpha _{n-2}}P_{n-2}(\xi ), \end{aligned}$$

(A.24a)

$$\begin{aligned} \alpha _{n-1}= & {} \frac{\beta }{n}\sqrt{\frac{4 n^2-1}{(\beta +1)^2 - n^2}}, \quad n=1,\ldots ,M \quad (\text {}M \le \beta ), \end{aligned}$$

(A.24b)

with $P_{-1}(\xi )=0$, $P_0(\xi )=\frac{1}{\sqrt{\beta +1}}$ and $\alpha _{-1}=1$ [5].

Using the Gram polynomials $P_n(\xi )$[5] with degree $M (M\le \beta )$, $f_{\text {match}}(z^{\text {left}})$ and $f_{\text {match}}(z^{\text {right}})$, we can obtain the projection coefficients such that

$$\begin{aligned} b_n^{\text {left}} = \sum _{i=0}^\beta f_{\text {match}}\left( z_i^{\text {left}}\right) P_n(\xi _i), \quad b_n^{\text {right}} = \sum _{i=0}^\beta f_{\text {match}}\left( z_i^{\text {right}}\right) P_n(\xi _i), \end{aligned}$$

(A.25)

where $z_i^{\text {left}}=b-\delta +i\Delta x$, $z_i^{\text {right}}=b+d+i\Delta x$ and $\xi _i=-1+\frac{2i\Delta x}{\delta }, i=0,\ldots ,\beta$.

To obtain values of the smooth matching function $f_{\text {match}}(z)$ in the extend domain $I_1=\left[ b,b+d\right]$ that satisfies the matching condition (13) and to implement efficiently, we use even and odd decompositions and obtain two smaller linear systems [15]. By splitting the original problem into smaller systems, it helps improve the conditioning of the original system. We note that a more flexible and efficient procedure of the discrete periodic extension method can be found in [2]. In this paper, however, we used even odd decompositions found in [15]. Define the even functions $f_n^{\text {even}}(z)$ and the odd functions $f_n^{\text {odd}}(z)$.

First, introduce a parameter $K(K<Q, \text {for an integer } Q \text { much larger than }\beta )$ to enhance accuracy in approximating the Gram polynomials and define the index set t(K) as

$$\begin{aligned} t(K) = {\left\{ \begin{array}{ll} \{i\in \mathbb {N} \mid -K/2+1 \le i \le K/2\}, \quad K \text { is even}\\ \{i\in \mathbb {N} \mid -(K-1)/2 \le i \le (K-1)/2\}, \quad K \text { is odd} \end{array}\right. }. \end{aligned}$$

(A.26)

Split this index set t(K) into an even index set $t_{\text {even}}(K)=\{j\in t(K) \mid mod(j,2)=0\}$ and an odd index set $t_{\text {odd}}(K)=\{j\in t(K) \mid mod(j,2)=1\}$. With these two index sets, we can define the even functions $f_n^{\text {even}}(z)$ and the odd functions $f_n^{\text {odd}}(z)$ on uniform grid points $z^{extend}_i = b+i\Delta x, \; i=1,\ldots ,\gamma -1$, such that,

$$\begin{aligned} f_n^{\text {even}}(z^{extend})=\sum _{k\in t_{\text {even}}(K)}\hat{a}_k^n e^{i \frac{\pi k}{\delta +d}z^{extend} },\quad f_n^{\text {odd}}(z^{extend})=\sum _{k\in t_{\text {odd}}(K)}\hat{b}_k^n e^{i \frac{\pi k}{\delta +d}z^{extend} }. \end{aligned}$$

(A.27)

Here the coefficients $\{ \hat{a}_k^n,\hat{b}_k^n \}$ are obtained by solving the overdetermined linear systems

$$\begin{aligned} {\left\{ \begin{array}{ll} \sum _{k\in t_{\text {even}}(K)}\hat{a}_k^n e^{i \frac{\pi k}{\delta +d}z_l } = P_n(\xi _l)\\ \sum _{k\in t_{\text {even}}(K)}\hat{a}_k^n e^{i \frac{\pi k}{\delta +d}z_r } = P_n(\xi _r) \end{array}\right. },\quad {\left\{ \begin{array}{ll} \sum _{k\in t_{\text {odd}}(K)}\hat{b}_k^n e^{i \frac{\pi k}{\delta +d}z_l }= P_n(\xi _l)\\ \sum _{k\in t_{\text {odd}}(K)}\hat{b}_k^n e^{i \frac{\pi k}{\delta +d}z_r }= -P_n(\xi _r) \end{array}\right. }, \end{aligned}$$

(A.28)

where $z_l = b-\delta +l\Delta z, \; l=0,\ldots ,Q$, $z_r = b+d+r\Delta z, \; r=0,\ldots ,Q$, $\xi _l = -1+ \frac{2 l \Delta z}{\delta }, \; l=0,\ldots ,Q, \; \xi _r = -1+ \frac{2 r \Delta z}{\delta }, \; r=0,\ldots ,Q$. Finally, the values of the matching function $f_{\text {match}}$ on $z^{extend} \in \left[ b,b+d\right]$ can be obtained as

$$\begin{aligned} f_{\text {match}}(z^{extend})= \sum _{n=0}^M \left[ \frac{b_n^{\text {left}} +b_n^{\text {right}} }{2}f_n^{\text {even}}(z^{extend}) + \frac{b_n^{\text {left}} -b_n^{\text {right}} }{2}f_n^{\text {odd}}(z^{extend})\right] . \end{aligned}$$

(A.29)