A Linear-complexity Tensor Butterfly Algorithm for Compressing High-dimensional Oscillatory Integral Operators

This paper presents a multilevel tensor compression algorithm called tensor butterfly algorithm for efficiently representing large-scale and high-dimensional oscillatory integral operators, including Green's functions for wave equations and integral transforms such as Radon transforms and Fourier transforms. The proposed algorithm leverages a tensor extension of the so-called complementary low-rank property of existing matrix butterfly algorithms. The algorithm partitions the discretized integral operator tensor into subtensors of multiple levels, and factorizes each subtensor at the middle level as a Tucker-type interpolative decomposition, whose factor matrices are formed in a multilevel fashion. For a $d$-dimensional integral operator discretized into a $2d$-mode tensor with $n^{2d}$ entries, the overall CPU time and memory requirement scale as $O(n^d)$, in stark contrast to the $O(n^d\log n)$ requirement of existing matrix algorithms such as matrix butterfly algorithm and fast Fourier transforms (FFT), where $n$ is the number of points per direction. When comparing with other tensor algorithms such as quantized tensor train (QTT), the proposed algorithm also shows superior CPU and memory performance for tensor contraction. Remarkably, the tensor butterfly algorithm can efficiently model high-frequency Green's function interactions between two unit cubes, each spanning 512 wavelengths per direction, which represents over $512\times$ larger problem sizes than existing butterfly algorithms. On the other hand, for a problem representing 64 wavelengths per direction, which is the largest size existing algorithms can handle, our tensor butterfly algorithm exhibits 200x speedups and $30\times$ memory reduction comparing with existing ones. Moreover, the tensor butterfly algorithm also permits $O(n^d)$-complexity FFTs and Radon transforms up to $d=6$ dimensions.