A fast butterfly algorithm for generalized Radon transforms

Low-rank approximations

Clearly the range and possibly other factors such as gradient of phase $\Phi(\mathbf{x},\mathbf{k})$ determine the degree of oscillations in the kernel $e^{ 2\pi i\Phi(\mathbf{x},\mathbf{k})}$ . Let $N$ be an integer power of two, which is on the order of the maximum value of $\vert\Phi(\mathbf{x},\mathbf{k})\vert$ for $\mathbf{x}\in X$ and $\mathbf{k}\in K$ (the exact choice of $N$ depends on the desired efficiency and accuracy of the algorithm, which will be made specific in numerical examples). The design of the fast algorithm relies on the key observation that the kernel $e^{ 2\pi i\Phi(\mathbf{x},\mathbf{k})}$ , when properly restricted to subsets of $X$ and $K$ , admits accurate and low-rank separated approximations. More precisely, if $A$ and $B$ are two square boxes in $X$ and $K$ , with sidelengths $w(A)$ , $w(B)$ obeying $w(A)w(B)\leq1/N$ -- in which case the pair $(A,B)$ is called admissible -- then

$$\left\vert e^{2\pi i \Phi(\mathbf{x},\mathbf{k})} -\sum_{t=1}^{r_... ...n} \alpha_t^{AB}(\mathbf{x})\beta_t^{AB}(\mathbf{k})\ \right\vert\leq \epsilon, \ \ \mathbf{x}\in A, \ \ \mathbf{k}\in B,$$ (8)

where $r_{\epsilon}$ is independent of $N$ for a fixed error $\epsilon$ . Here and below the subscript $t$ is slightly abused: $t$ should be understood as multi-indices $(t_1,t_2)$ , and accordingly $r_{\epsilon}$ is the total number of terms in a double sum. Furthermore, Candès et al. (2009) showed that this low-rank approximation can be constructed via a tensor-product Chebyshev interpolation of $e^{ 2\pi i\Phi(\mathbf{x},\mathbf{k})}$ in the $\mathbf {x}$ variable when $w(A)\leq 1/\sqrt{N}$ , and in the $\mathbf {k}$ variable when $w(B)\leq 1/\sqrt{N}$ .

Specifically, when $w(B)\leq 1/\sqrt{N}$ , $\alpha_t^{AB}$ and $\beta_t^{AB}$ are given by

$$\alpha_t^{AB}(\mathbf{x})=e^{2\pi i \Phi(\mathbf{x},\mathbf{k}_t^B)},$$ (9)

$$\beta_t^{AB}(\mathbf{k})=e^{-2\pi i \Phi(\mathbf{x_0}(A),\mathbf{k}_t^B)}L_t^B(\mathbf{k}) e^{2\pi i \Phi(\mathbf{x_0}(A),\mathbf{k})};$$ (10)

and when $w(A)\leq 1/\sqrt{N}$ , $\alpha_t^{AB}$ and $\beta_t^{AB}$ are given by

$$\alpha_t^{AB}(\mathbf{x})=e^{2\pi i \Phi(\mathbf{x},\mathbf{k}_0(B))}L_t^A(\mathbf{x}) e^{-2\pi i \Phi(\mathbf{x}_t^A,\mathbf{k}_0(B))},$$ (11)

$$\beta_t^{AB}(\mathbf{k})=e^{2\pi i \Phi(\mathbf{x}_t^A,\mathbf{k})}.$$ (12)

Boldface letters $\mathbf{k}_t^B$ , $\mathbf{x}_t^A, \mathbf{k}_0(B), \mathbf{x}_0(A)$ denote 2D vectors. $\mathbf{k}_t^B$ is a point on the 2D, $q_{k_1}\times q_{k_2}$ Chebyshev grid in box $B$ centered at $\mathbf {k}_0(B)$ , i.e., let $\mathbf{k}_t^B=(k_{t_1}^B,k_{t_2}^B)$ , $\mathbf{k}_0(B)=(k_{0_1}^B,k_{0_2}^B)$ , then

$$k_{t_1}^B=k_{0_1}^B+w(B)z_{t_1}, \quad 0\leq t_1 \leq q_{k_1}-1,$$ (13)

$$k_{t_2}^B=k_{0_2}^B+w(B)z_{t_2}, \quad 0\leq t_2 \leq q_{k_2}-1,$$ (14)

where

$$\left\{z_{t_i}=\frac{1}{2}\cos\left( \frac{\pi t_i}{q_{k_i}-1}\right)\right\}_{0\leq t_i\leq q_{k_i}-1,\ i=1,2}$$ (15)

is the 1D Chebyshev grid of order $q_{k_i}$ on $[-1/2,1/2]$ . See Figure 1 for an illustration. $L_t^B(\mathbf{k})$ is the 2D Lagrange interpolation defined on the Chebyshev grid:

$$L_t^B(\mathbf{k})=\left(\prod_{s_1=0, s_1\neq t_1}^{ q_{k_1}-1}\f... ...=0, s_2\neq t_2}^{ q_{k_2}-1}\frac{k_2-k_{s_2}^B}{k_{t_2}^B-k_{s_2}^B} \right).$$ (16)

Analogously, $\mathbf{x}_t^A$ is a point on the 2D, $q_{x_1}\times q_{x_2}$ Chebyshev grid in box $A$ centered at $\mathbf{x}_0(A)$ , and $L_t^A(\mathbf{x})$ is the 2D Lagrange interpolation defined on this grid. Based on the discussion above, the number $r_{\epsilon}$ in low-rank approximation 8 is equal to $q_{k_1}q_{k_2}$ when $w(B)\leq 1/\sqrt{N}$ , and $q_{x_1}q_{x_2}$ when $w(A)\leq 1/\sqrt{N}$ .

grid
Figure 1. A 2D, $q_{k_1}\times q_{k_2}$ ($q_{k_1}=7$ , $q_{k_2}=5$ ) Chebyshev grid in box $B$ . $\mathbf {k}_0(B)$ is the center of the box. $\mathbf {k}_t^B=(k_{t_1}^B,k_{t_2}^B), \ 0\leq t_1 \leq q_{k_1}-1, \ 0\leq t_2 \leq q_{k_2}-1$ is a point on the grid.

A simple way of viewing expressions 9 - 12 is: when $w(B)\leq 1/\sqrt{N}$ , plugging expression 9 into approximation 8 (leaving $\beta_t^{AB}(\mathbf{k})$ as it is) yields

$$e^{2\pi i \Phi(\mathbf{x},\mathbf{k})} \approx \sum_t e^{2\pi i \Phi(\mathbf{x},\mathbf{k}_t^B)}\beta_t^{AB}(\mathbf{k}), \ \ \mathbf{x}\in A, \ \ \mathbf{k}\in B.$$ (17)

For fixed $\mathbf {x}$ , the right hand side of equation 17 is just a special interpolation of function $e^{ 2\pi i\Phi(\mathbf{x},\mathbf{k})}$ in variable $\mathbf {k}$ , where $\mathbf{k}_t^B$ are the interpolation points, $\beta_t^{AB}(\mathbf{k})$ are the basis functions. Likewise, when $w(A)\leq 1/\sqrt{N}$ , plugging expression 12 into approximation 8, we get

$$e^{2\pi i \Phi(\mathbf{x},\mathbf{k})} \approx \sum_t e^{2\pi i \Phi(\mathbf{x}_t^A,\mathbf{k})}\alpha_t^{AB}(\mathbf{x}), \ \ \mathbf{x}\in A, \ \ \mathbf{k}\in B.$$ (18)

For fixed $\mathbf {k}$ , the right hand side of equation 18 is a special interpolation of $e^{ 2\pi i\Phi(\mathbf{x},\mathbf{k})}$ in variable $\mathbf {x}$ : $\mathbf{x}_t^A$ are the interpolation points, $\alpha_t^{AB}(\mathbf{x})$ are the basis functions.

Once the low-rank approximation 8 is known, computing the partial sum

$$u^{B}(\mathbf{x}):=\sum_{\mathbf{k}\in B}e^{2\pi i \Phi(\mathbf{x},\mathbf{k})}g(\mathbf{k})$$ (19)

generated by points $\mathbf {k}$ inside a box $B$ becomes

$$u^{B}(\mathbf{x})\approx \sum_{\mathbf{k}\in B} \sum_t \alpha_t^{... ...{AB}(\mathbf{k}) g(\mathbf{k}) = \sum_t \alpha_t^{AB}(\mathbf{x})\delta_t^{AB}, \ \ \mathbf{x}\in A,$$ (20)

where

$$\delta_t^{AB}:=\sum_{\mathbf{k}\in B} \beta_t^{AB}(\mathbf{k})g(\mathbf{k}).$$ (21)

The case that the box $B$ represents the whole domain $K$ is of particular interest, since it corresponds to the original problem. Therefore, if we can find the set of interaction coefficients $\delta_t^{AB}$ relative to all admissible couples of boxes $(A,B)$ with $B=K$ , our problem will be solved.

A fast butterfly algorithm for generalized Radon transforms