  Image subtraction is essential for transient detection in time-domain
astronomy. The point spread function (PSF), photometric scaling, and sky
background generally vary with time and across the field-of-view for imaging
data taken with ground-based optical telescopes. Image subtraction algorithms
need to match these variations for the detection of flux variability. An
algorithm that can be fully parallelized is highly desirable for future
time-domain surveys. Here we show the Saccadic Fast Fourier Transform (SFFT)
algorithm for image differencing. SFFT uses $\delta$-function basis for kernel
decomposition, and the image subtraction is performed in Fourier Space. This
brings about a remarkable improvement of computational performance of about an
order of magnitude compared to other published image subtraction codes. SFFT
can accommodate the spatial variations in wide-field imaging data, including
PSF, photometric scaling, and sky background. However, the flexibility of the
$\delta$-function basis may also make it more prone to overfitting. The
algorithm has been tested extensively in real astronomical data taken by a
variety of telescopes. Moreover, the SFFT code allows for the spatial
variations of the PSF and sky background to be fitted by spline functions.
