Partial Fourier transform of large sparse 3D array

I am trying to calculate the Fourier transform of a large 3D array (30,000 x 30,000 x 1500 pixels) with about 30,000 x 15,000 nonzero elements. There is at most 1 nonzero pixel in each x-y plane (first 2 dimensions). I cannot downsample the array because I need enough input resolution to get a specific data span in the output.

To avoid storing/loading multiple 10TB arrays, I'm thinking of doing 1D Fourier transform on each dimensions sequentially. I would appreciate further suggestions to speed up computation and reduce memory load. Moreover, if I only want to compute the first 256 x 256 x 768 pixels of the output, is there a way to further reduce computation?

[David Goodmanson]

Hi Linh,

you say that there is at most 1 nonzero pixel in each x-y plane, but there are only 1500 x-y planes so this appears to conflict with the statement that there are 30,000 x 15,000 nonzero elements.

Answer 1

[Matt J]

You can do the partial Fourier transforms by computing a set of partial DFT matrices:

dftPartial=@(N,k)sparse(exp(-(2j*pi/N)*(0:N-1).*(0:k-1)'));

Fxy=dftPartial(3e4,256);
Fz= dftPartial(1500,768);

and combine them into a 3D operator by downloading the KronProd class,

https://www.mathworks.com/matlabcentral/fileexchange/25969-efficient-object-oriented-kronecker-product-manipulation

F=KronProd({Fxy,Fz},[1,1,2])

The 3D sparse array, X, that is being transformed can be represented as an ndSparse array,

https://www.mathworks.com/matlabcentral/fileexchange/29832-n-dimensional-sparse-arrays?s_tid=srchtitle

Then, you shoul be able to do simply,

Y=full(F*X);

[Matt J]

This example took about 2 min. on my 16 GB laptop

>> X=ndSparse.sprand(F.domainsizes,1500/size(F,2));

>> tic; Y=full(F*X); toc

Elapsed time is 126.221071 seconds.

>> whos F X Y

Name Size Kilobytes Class Attributes

F 50331648x1350000000000 207247 KronProd

X 30000x30000x1500 36 ndSparse sparse

Y 256x256x768 786432 double complex

Answer 2

[Matt J]

There is at most 1 nonzero pixel in each x-y plane (first 2 dimensions) ... I'm thinking of doing 1D Fourier transform on each dimensions sequentially.

That would seem to be more work than necessary. The 2D DFT of an impulse (and any sub-rectangle of it) is something that can be computed anaytically. So, if each slice is an impulse image, it's really only the FFT along the 3rd dimension that you would need to do numerically, which should be very tractable.