So I'm following this example on how to compress a picture with compressed sensing.

http://www.pyrunner.com/weblog/2016/05/26/compressed-sensing-python/

In the example, the picture is sparse in the Frequency domain, so the example implements a matrix that includes a DCT. That Matrix is in charge of the DCT and includes the sensing matrix.

So if we typically imply: $y=\phi x$ where $\phi$ is the sensing matrix, $x$ the samples and $y$ the resulting compressed signal.

He implies: $b=Ax$ and $A\equiv \phi \psi$ and $\psi$ converts from the spectral domain to the original domain. $b$ should, therefore, be the output ($y$).

Ok, I'm trying to apply this example to a fake case where the image is sparse in the original domain.

So I don't really need to convert to the frequency domain. Therefore I guess that my fake example shall be easier.

So I eliminated the DCT part, the Kronecker product and so on (in the example it is detailed that part). But then I got quite lost.

Maybe it's a problem of not properly understanding the process from the theoretical point of view. So, It is my understanding that, in order to apply the compressed sending, the process is as follows:

  1. Find the proper $\phi$ (or in my case $A$ without the $\psi$ so $A=\phi$). Ensuring the necessary properties (eg: check the coherence of the sensing matrix $\mu(\phi)$, at the end of the day you have to calculate the matrix and ensure that it will not map two inputs ($x_1$ and $x_2$, where $x_1\ne x_2$) to the same $y$).
  2. Convert your input, with the fresh new calculated sensing matrix
  3. Recover the original signal from $y$. For instance finding the signal $x$ that minimizes the $\ell_1$ norm.

Now, in the example, in a particular moment he creates the matrix $A$, as follows:

# create dct matrix operator using kron (memory errors for >large ny*nx)
A = np.kron(
    spfft.idct(np.identity(nx), norm='ortho', axis=0),
   spfft.idct(np.identity(ny), norm='ortho', axis=0)
   )
A = A[ri,:] # same as phi times kron

As I don't really need all the DCT part I just create an identity matrix nx*ny, and then I sample it with ri which is basically a group of randomly selected indexes (this is how the sampling process of the original picture is done).
Cool, first problem: he makes an IDCT of the identity matrix and then he has values in the matrix, but as I pick random samples from identity matrix directly, most of the original sensing matrix is 0 values.

So:

  1. Most of my sensing matrix is 0
  2. I did not check that the matrix is not mapping two different inputs to the same value. Indeed if the Matrix is almost 0 (or all 0), it will. Indeed it will "kill" whatever I use as an input $x$, as $y$ will always be 0.

On top of that, I don't see how it calculates the sensing matrix.

He does:

# do L1 optimization
vx = cvx.Variable(nx * ny)
objective = cvx.Minimize(cvx.norm(vx, 1))
constraints = [A*vx == b]
prob = cvx.Problem(objective, constraints)
result = prob.solve(verbose=True)
Xat2 = np.array(vx.value).squeeze()

Is he using as the variable vx? But vx is $x$, right? so, at the end of the day, is he trying to find the input that, once multiplied by the sensing matrix provides you with the output? Does it make sense?

I guess I'm misunderstanding the meaning of the above code and I don't really get what cvx is doing there. Because I don't get how this will help to set A and the calculate the compressed output.

So, to summarize:

  1. Is the general process, from the theoretical point of view, in general as described above?
  2. Am I wrong when I say that he is trying to find using cvx the input that once multiplied by A provides you with the output? I guess he has to calculate the A matrix that will not map two different X input to the same b output, then calculate the b value with the sensing matrix and the provided input.
  • Try with a randomly generated sensing matrix, instead of designing it yourself. To my understanding your sampling matrix is also sparse which is problematic when you are not using DCT. DCT of image actually spreads the non zero coefficients. – MimSaad Oct 15 '17 at 20:43
  • Hi, thanks a lot for your help. So, I've been checking the Sensing matrix out. As you point out, it is sparse, as it is mostly 0. I noticed that I needed to increase the number of samples that compose the "ri" index in order to obtain good results (it also depends on the signal of interest degrees of freedom, right?). I checked it out and from an original sensing matrix of 6724x6724 I end up with a matrix of 6052x6724 (90%). It works with a little bit lower value. – f.gallardo Oct 16 '17 at 10:26
  • I checked it out and I have enough number of ones in the sensing matrix at that point. So that is why it is working (maybe), and I think it is fine in terms of reconstruction as the spark should be high enough as I'm using the identity matrix, right? Now, that means that maybe I don't really need to find the Sensing matrix, as I was stating before. That is good. I was misunderstanding that part of the process. – f.gallardo Oct 16 '17 at 10:26
  • On top of that, after reviewing the code I noticed that "b" is the sampled input, not the output. So, when he makes: "constraints = [A*vx == b] " he is trying to find a "vx" that, taking into account the used sensing matrix, makes the product of A (the sensing matrix) times "vx" as similar as possible to b (the sampled input). – f.gallardo Oct 16 '17 at 10:31
  • In principle every and each element of Y vector (measurements) should contain some information about the most of signal elements. When you are using an sparse sampling matrix for an sparse signal (a signal that is sparse in the domain that the sensing matrix is sparse) directly, you miss a lot of non-sparse coefficients by multiplying them in zero, that is why when your increased number of 1's you achieved better results. – MimSaad Oct 16 '17 at 19:36

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.