6
$\begingroup$

I have been posting this kind of question in Cross Validated, but since this one deals almost entirely with mathematics, I will post it here.

In signal reconstruction using compressed sensing, we want to sample a signal $f$ in order to obtain a smaller (compressed) signal $b$:

$$b = \phi f$$

However, $f$ can also be expressed by a linear combination of basis functions $\Psi$ and its coefficients $c$:

$$f = \Psi c$$

So the first equation and second equation together produce:

$$b = \phi \Psi c$$

Furthermore, we also want to promote the sparsity of $c$. This is done by solving the following optimization problem:

$$\text{min } ||c||_{l_{1}} \text{ subject to } b = \phi \Psi c$$

where $||c||_{l_{1}}$ is the $l_{1}$-norm.

I think I understand this part but I'd like to know more about the bases $\Psi$. How are they chosen? How are they implemented in an algorithm? Right now I'm trying to understand a couple of examples that use the following $\Psi$ and $\Phi$:

1.

D=dct(eye(n,n)); % \Psi
A=D(perm,:); % \Phi \Psi

where dct is the discrete cosine transform, n is the dimension of the original signal and perm are the first numbers of a list of randomly generated numbers:

r1 = permutation(arange(1, n))
perm = r1[0:m]

2.

Atest = zeros((nx, ny)).reshape(1, nx*ny)
Adelta = zeros((k, nx*ny))
for i, j in enumerate(r1k):
    Atest[0,j] = 1
    Adelta[i, :] = dct(Atest)
    Atest[0, j] = 0

here nx and ny are the dimensions of the original signal, r1k is also permutation (similar to perm). Adelta is produced by choosing a point in a matrix Atest, transform it using dct and adding the result as a row in matrix Atest.

In these pieces of code, I know A and Adelta represent $\Phi \Psi$ but I don't really understand why. I'd appreciate a few comments about this.

By the way, if you want to take a look at the full example, here is one using MATLAB and this one in Python. This example is also related to my other question in Cross Validated, so if you want to contribute to that also, you are more than welcome, although it deals with an implementation issue which may not be interesting enough.

UPDATE:

An additional question: How should we deal with the actual reconstruction. According to the second equations $f = \Psi c$, once we know $c$, an approximation of $f$ is found simply by calculating $\Psi c$. However, what I see in practice is actually the direct application of a DCT to the coefficients $c$. Why?

Thanks

$\endgroup$
  • $\begingroup$ For part 1 above, $D$ is $\Psi$ and $A$ is $\Phi \Psi$. The value of $\Phi$ is implied by the selection of the random rows of $D$, so it's really just a row-selection matrix. For instance, if perm(1)=13, then the first row of $\Phi$ is all zeros except with a 1 in column 13. Likewise, if perm(2)=47, then the second row of $\Phi$ is all zeros except a 1 in column 47. Does this answer your question? $\endgroup$ – AnonSubmitter85 Jun 23 '13 at 0:51
  • $\begingroup$ Actually, it's the rows that you're selecting based on perm. I think you're also missing the dct function that is being applied on the rows before such selection takes place. $\endgroup$ – Robert Smith Jun 23 '13 at 2:35
  • $\begingroup$ I understand that perm is selecting rows, that's why it implies that $\Phi$ is a row-selection matrix. As for the DCT function, taking the DCT of an identity matrix is just an east way to create the DCT matrix. The following are equivalent: dct(X) == dct(eye(size(X))) * X. $\endgroup$ – AnonSubmitter85 Jun 23 '13 at 2:40
  • $\begingroup$ Are you confused about what A represents? If so, in the first case, it simply represents the values of the DCT that you have available, that is, your incomplete information from which you are going to reconstruct the entire signal. The DCT is a linear operator and as such has a matrix representation, thus, the line D = dct(eye(n,n)) is simply creating the matrix representation of the DCT. Removing rows from D is the same thing as removing values from the DCT of the signal. Thus, A * x is equal to values of dct(x) at the indices in perm. $\endgroup$ – AnonSubmitter85 Jun 23 '13 at 2:50
  • $\begingroup$ I don't think I'm confused about $A$. It is mostly "How are $\Psi$ and $\Phi$ chosen?" and "How are they implemented in an algorithm?" The second one might not be a good fit for this site, though. $\endgroup$ – Robert Smith Jun 23 '13 at 3:25
2
$\begingroup$

I will not enter in the details of your code (which is not as readable as it should be), i will rather give you an abstract description: $\phi$ is the matrix that represents the way you sample your signal $f$. For example, if you are doing subsampling, i.e. you just sample a subset of the entries of $f$, then $\phi$ is the identity matrix with rows removed. On the other hand $\Psi$ represents the basis that you choose to expand your signal $f$. For example, if you know that $f$ is a time series containing only a couple of frequencies, then you can take $\Psi$ to be the Discrete Fourier (or discrete cosine) matrix and then you try to solve for $c$. Note that $c$ will be sparse, with nonzero entries associated to the distinct frequencies.

$\endgroup$
  • $\begingroup$ Thanks. Let me know what I can do to make the code readable (I didn't want to post the whole code because it would detract from the main question.) In any case, yes, that was my understanding of $\Psi$ and $\Phi$. Maybe the issue here is that I don't understand what leads to choose an operation such as DCT(I(n, n)) unlike the example you described in which it is very clear the purpose. $\endgroup$ – Robert Smith Jun 22 '13 at 23:41
  • $\begingroup$ I would assume that this is because the signal that you the designer of this example wants to sample is sparse in frequency. But which one is this signal $f$ in the code you show? I recommend this, and its quite easy: why don't you design an example in MATLAB exactly along these lines? I.e. construct a time series with e.g. two frequencies, take some samples (b vector) and then using the DCT matrix solve the l1 problem and see if you can recover your time series. $\endgroup$ – Manos Jun 22 '13 at 23:50
  • $\begingroup$ I already did that (in Python) from an example in MATLAB (mathworks.com/company/newsletters/articles/…) It works but I still find rather mysterious the basis part. Actually, most of what I have done so far can be found in this other question (stats.stackexchange.com/questions/62280/…) $\endgroup$ – Robert Smith Jun 22 '13 at 23:56
  • $\begingroup$ By the way, how is your linear algebra? $\endgroup$ – Manos Jun 23 '13 at 1:11
  • 1
    $\begingroup$ It all depends on what your signal is. That's the key to what basis one has to use. For example, in image analysis, its very common to use $\Psi$ as the wavelet transform. Why? Because usually images are sparse in the wavelet domain. So you should really think what is your signal $f$ in your example or your application. $\endgroup$ – Manos Jun 23 '13 at 2:31
0
$\begingroup$

To take case 1, you know that $A=\Phi\Psi$. The DCT matrix, $D$, is $\Psi$. The value of $\Phi$ is implied by the random selection of rows in the line A = D(perm,:);. This is equivalent to multiplying D by a row-selection matrix: A = P * D. To give a specific example, say that D is a 5x5 matrix and you want to make the 2x5 matrix A that is made up of the rows 4 and 2 of D. Then

$$ P= \left[ \begin{array}{ccccc} 0 & 0 & 0 & 1 & 0\\ 0 & 1 & 0 & 0 & 0\\ \end{array} \right] $$

The following are good background material on the theoretical aspects of the problem:

Uncertainty Principles and Signal Recovery

Optimally Sparse Representation in General Nonorthogonal Dictionaries

Sparsity and Incoherence in Compressive Sampling

Decoding by Linear Programming

Uncertainty Principles and Ideal Atomic Decomposition

$\endgroup$
  • $\begingroup$ @RobertSmith To the contrary, they cover it in exacting detail. Also, I would not recommend reading a single part of these papers in isolation. $\endgroup$ – AnonSubmitter85 Jun 23 '13 at 0:04
  • $\begingroup$ Then would you mind to point out the page in which it's covered in exacting detail? $\endgroup$ – Robert Smith Jun 23 '13 at 0:05
  • $\begingroup$ Pages 2-17 of "Optimally Sparse Representation in General Nonorthogonal Dictionaries". $\endgroup$ – AnonSubmitter85 Jun 23 '13 at 0:08
  • $\begingroup$ @RobertSmith You do realize that dictionary == bases, right? $\endgroup$ – AnonSubmitter85 Jun 23 '13 at 0:10
  • $\begingroup$ Right. However, there are only a couple of paragraphs that are directly applicable to my question. $\endgroup$ – Robert Smith Jun 23 '13 at 0:13
0
$\begingroup$

You say "but I'd like to know more about the bases $\Psi$. How are they chosen? How are they implemented in an algorithm? "

$\Psi$ is chosen so that $c$ is sparse. If $\Psi$ is equal to the identity, then the signal is sparse in the canonical basis.

Let us imagine the signal is sparse in the canonical basis (aka $\Psi =\mathbf{I}$) then the difference between a simple sampling problem and compressive sensing problem is the difference between $\Phi$ being the identity and an acceptable compressive sensing measurement matrix.

How is it implemented in algorithm? $\Psi$ is either a matrix made up of columns of the discretized version of each function of the basis or just a subroutine.

$\endgroup$

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.