A fast butterfly algorithm for generalized Radon transforms

Numerical complexity and accuracy

To analyze the algorithm's numerical complexity, let us assume the numbers of Chebyshev points in every box and every dimension of $K$ and $X$ are all equal to a small constant $q$ , i.e., $q_{k_1}=q_{k_2}=q_{x_1}=q_{x_2}=q$ and $r_{\epsilon}\equiv q^2$ . The main workload of the fast butterfly algorithm is in Steps 2 and 4. For each level, there are $N^2$ pairs of boxes $(A,B)$ , and the operations between each $A$ and $B$ is $O(r_{\epsilon}^2)$ , which can be further reduced to $O(r_{\epsilon}^{3/2})$ by performing Chebyshev interpolation one dimension at a time. Since there are $\log N$ levels, the total cost is $O(r_{\epsilon}^{3/2} N^2\log N)$ . It is not difficult to see that Step 3 takes $O(r_{\epsilon}^2 N^2)$ , and Steps 1 and 5 take $O(r_{\epsilon}N_fN_h)$ and $O(r_{\epsilon}N_{\tau}N_p)$ operations. Considering the initial Fourier transform of preparing data in the $(f,h)$ domain, we conclude that the overall complexity of the algorithm is $O(N_h N_t \log N_t+ r_{\epsilon}^{3/2}N^2\log N +r_{\epsilon}^2 N^2+ r_{\epsilon}(N_fN_h+N_{\tau}N_p))$ . The analysis in Candès et al. (2009) shows that the relation between $r_{\epsilon}$ and error $\epsilon$ is $r_{\epsilon}=O(\log ^4(1/\epsilon))$ . We would like to mention that this is only the worst case estimate. Numerical results in the same paper demonstrate that the dependence of $r_{\epsilon}$ on $\log(1/\epsilon)$ is rather moderate in practice.

In comparison, the conventional velocity scan requires at least $O(N_{\tau}N_pN_h)$ computations, which quickly becomes a burden as the problem size increases. Yet the efficiency of our algorithm is mainly controlled by $O(N^2\log N)$ with a constant polylogarithmic in $\epsilon$ , where $N$ depends neither on data size nor on data content (here we mean the data after the Fourier transform). Since the Chebyshev interpolation is only performed on the kernel, our choice of parameters ($N$ and number of Chebyshev points) relies on the preknowledge about the range of $f$ , $h$ , $\tau$ , and $p$ . In other words, we need a general idea about how oscillate the kernel is. Recall that everything is mapped to a unit square, so the larger the range of $\Phi(\mathbf{x},\mathbf{k})$ is, the more oscillations occur in the unit square. If the original data (data before the Fourier transform) contain high frequency information, the accuracy will be affected as the frequency bandwidth is now larger. A possible way to get around it is to divide the Fourier domain into two or three smaller subdomains (so the range of $f$ in each subdomain is smaller than the original problem), and apply the fast algorithm to each part separately, finally add the results back together. This only increases the cost by a small factor, but presumably offers better accuracy.

A fast butterfly algorithm for generalized Radon transforms