A Superfast Direct Solver for 2D Type-II Inverse Nonuniform Discrete Fourier Transform Based on Hierarchically Semiseparable Matrices

Abstract

Reconstructing images on Cartesian grids from nonuniform Fourier samples is a basic computational task in Fourier reconstruction methods for computed tomography. This paper proposes a direct solver for the 2D type-II nonuniform discrete Fourier transform~(NUDFT) based on the face-splitting product of hierarchically semiseparable~(HSS) matrices. We define the face-splitting product of two HSS trees and prove that the face-splitting product of two HSS matrices admits an HSS representation on the resulting tree. The proof is constructive, providing explicit formulas for the HSS generators and bounds on the off-diagonal ranks at each level. For the NUDFT, a Fourier change of basis expresses the transformed 2D matrix as a face-splitting product of two transformed 1D NUDFT matrices. This identity yields an algebraic construction of a 2D HSS approximation from two 1D HSS approximations, for which we establish a Frobenius-norm approximation error bound. Combining this construction with recompression, URV factorization, and a two-dimensional inverse fast Fourier transform gives a direct solver for the associated least-squares problem. For fixed accuracy, square frequency grids, and balanced HSS trees with fixed leaf column sizes, the solver has an offline complexity of $O(M{\log }^{2}N+N^{3/2}{\log }^{3}N)$ and an online complexity of $O(M+N{\log }^{3}N)$ per right-hand side, where $M$ and $N$ are the numbers of sample points and unknown coefficients, respectively. Numerical experiments on random and polar grids illustrate the computational scaling and reconstruction accuracy of the method, and demonstrate a substantial reduction in lsqr iteration counts when the direct solver is used as a preconditioner.