3
$\begingroup$

When I use the cvx matlab toolbox, I met a puzzled problem. The function of fft (or dct, wavelet, etc.) cannot be recognized by the type of 'cvx'. For the 1-d fft, it can be constructed to an equivalent matrix. But for the 2-d fft, how to transform it to the matrix form ?

$\endgroup$
5
$\begingroup$

Yes, it is true, the CVX package (disclosure: I am the author) does not support the use of the fft() command in constraint and objective expressions. However, it's relatively simple to get the equivalent result, in both the 1D and 2D cases.

First, 1D. Given a vector $x\in\mathbb{R}^n$, the DFT of $x$ is equal to $Wx$, where $W$ is the so-called DFT matrix. Here are a few ways to construct the DFT matrix in MATLAB:

W = fft(eye(n))
W = exp(1j*2*pi*[0:n-1]'*[0:n-1]/n)
W = exp(bsxfun(@times,1j*2*pi*[0:n-1]',[0:n-1]/8))

Now try this with random vectors to verify correctness:

x=randn(n,1)+1j*randn(n,1)
fft(x)
W*x

What about 2D ffts? Given a square matrix $X\in\mathbb{R}^{n\times n}$, the 2D DFT can be obtained by applying the 1-D DFT to each column of the matrix, and then applying the 1-D DFT to each row of the matrix. Expressed in terms of the DFT matrix $W$, the 2D DFT is simply $WXW^T$. Go ahead, try the following MATLAB expressions and see that they are equivalent:

X = randn(n,n)+1j*randn(n,n)
fft2(X)
fft(fft(X,[],2),[],1)
fft(fft(X).').'
W*X*W.'

Note the use of the non-conjugated transpose W.' above; that's important. If you use the standard Hermitian transpose instead, ', you'll get a scrambled result.

In summary, given a square matrix variable X, you can perform a 2D FFT on a square CVX variable as follows:

W=fft(eye(size(X));
W*X*W.'

Again, note that I am using a non-conjugated transpose above.

Please note that 2D DFTs get big fast. That is, you have n^2 variables here if the dimensions of your matrix are [n n]. So even a 512 by 512 image, not that large by matrix standards, is still over 250k variables, and that can be a lot by CVX standards. Once you are satisfied with your model's performance on smaller problems, you may need to develop your own numerical engine to solve larger cases.

Given that it is relatively simple to work around CVX's limitation here, it is reasonable to ask: why not just add support for fft() in CVX, using one of techniques described here? Well, I may do that. However, it's important to remember that FFT stands for fast Fourier transform. It is a fast algorithm for generating the discrete Fourier transform (DFT). Specifically, it can compute the DFT in only $O(N\log N)$ operations.

The matrix workaround, on the other hand, constructs a full DFT matrix at a cost of $O(N^2)$ operations---so much for "fast"! I do not want to deceive users with a false sense that CVX is really using the fast Fourier transform. The fact is, it cannot: the underlying solvers do not exploit the $N\log N$ efficiency of the FFT operator.

The CVX example library contains a number of problems that utilize a manually generated Fourier transform. In particular, look at the filter design examples.

$\endgroup$
  • $\begingroup$ Thanks Michael! How about for 2-d wavelet or curvelet ? It seems like not so easy to be transformed to some matrices. $\endgroup$ – syli Apr 24 '13 at 1:52
  • $\begingroup$ I am afraid I do not know. You will have to consult the math for each transform separately. But my suspicion is than once you form the 1D matrix, the 2D result is similar to this. $\endgroup$ – Michael Grant Apr 24 '13 at 2:34
1
$\begingroup$

Say A is a M by N matrix, then fft2(A) = dftmtx(M) * A * dftmtx(N). Also, we have $$\text{vec}(ABC)=(C^T \otimes A)\text{vec}(B)$$ So you can use Kronecker product to construct the 2D DFT matrix. See here for details.

$\endgroup$
  • $\begingroup$ I'm glad you turned your comment into an answer---I voted it up. With CVX, the Kronecker approach is not necessary, because you can work with variables in their matrix form. But in other contexts this is a crucial step. In effect, this is what CVX has to do internally. $\endgroup$ – Michael Grant Apr 23 '13 at 19:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.