# Higher-order Kronecker product

I am trying to generate a 2D DFT matrix in matlab, which I need for 2D compressed sensing (CS) problems.

Lets say $N_1=8$, $N_2=16$, then according to the requirement of CS, to generate a 2D DFT matrix, we need to first generate DFT matrices $A$ $(N_1\times N_1)$ and $B$ $(N_2\times N_2)$. Later calculate Kronecker product of $A$ and $B$ which will be of order $N_1\times N_2\times N_1 \times N_2$.

A= dftmtx(N1);
B= dftmtx(N2);dft2D = kron(A,B);


It can be done by Matlab, but there is a limit to the max size one can calculate (in my case its $N_1=64$, $N_2=256$). But I need it for higher orders i.e., $N_1=128$, $N_2=512$ or higher orders. Is there any efficient way to do this without memory issues?

• Which solver are you using? If you use a solver that supports supplying the matrix (and probably its adjoint) as functions instead of as matrices, you can implement the operators via separate FFTs along each dimension. This will result in huge savings in computational complexity. Apr 12 '18 at 8:06
• Do you have a preferred language for programming this - Python, Matlab? Apr 12 '18 at 8:38
• @thomas arildsen I am using matlab only. Apr 12 '18 at 9:01

Since you are using Matlab, I suggest using SPGL1 as your reconstruction solver. It can handle complex numbers and you can specify your matrix operator as a function.

This example demonstrates how to specify your matrix operator as a function instead of a matrix (see opA).

In your case, I would define your operator as a function that:

1. takes a vector x as input,
2. reshapes that vector to a matrix X of your desired shape,
3. applies fft to the matrix: Z = fft(X) (Matlab's fft transforms each column separately if you give it a matrix),
4. applies fft to the transpose of the result: W = fft(Z'),
5. reshapes W back into a vector w and returns that

If your operation additionally includes some sub-sampling (which it usually does in order to be compressed sensing) then that can be implemented as indexing into w for the entries that are sampled.

The above steps imply that x is in the domain of the solution you are seeking and therefore sparse. Thus, the sparse vector here is in the time domain and you are sampling in the Fourier domain. If, on the other hand you want to model a signal that is sparse in frequency and you sample in the time domain, you should replace fft by ifft above.

• I am not so familiar with matlab operator as function coding. I tried understanding the spgl1 solver but not able to learn it. I don't have much time remaining to finish my thesis. Can you send me the code. for an image of 1024*256 dimension, image being already in frequency domain. and use 700 samples in each column whose index is predefined for each column. May 3 '18 at 14:16
• No, I am not interested in spending my time coding it for you. I have tried to help as well as I reasonably could. It ought to be possible to implement it from the description above without too much trouble. May 9 '18 at 20:18

You want to wrap the call to fft in a function that does the reshaping and whatnot for you. Under only the most trivial of circumstances do you ever want to explicitly construct the matrix. With the 2-D DFT as an example, a 1024 x 1024 input array would require a 1024^2 x 1024^2 matrix. Moreover, the matrix implementation will be an order of magnitude slower than the FFT implementation.

For instance, suppose you are using SPGL1, where the format of the function handle is f(x,mode), where if mode == 1, the forward linear operator is applied to the input, and if mode == 2, the adjoint is applied.

function y = fft2_operator(x,mode,M,N)
x = reshape(x,[M,N]);  % Reshape the input to an M x N array.
if(mode == 1)
y = fft2(x);
else
y = ifft2(x);
end
y = y(:);  % We want to return an MN x 1 vector.
end


You can then define a function handle to pass to the solver using something like the following:

% The size of the input with some made up data.
M = 1000;
N = 1024;
x = rand(M,N);

f_fft = @(x,mode) fft_operator(x,mode,M,N);
y = spgl1(f_fft,x);


If you are using a different solver, then the details of the above will change but the approach will be the same.