Compressed Sensing : A python demo

[Dilawar Singh]

TL;DR. Dont read. Only useful for signal/image-processing/information theory people.

Can you solve Ax = b for x when number of rows are less than number of columns in A i.e there are more variables than equations? Such a system is called under-determined system. For examples, give me 3 numbers whose average is 3. Or solve x/3 + y/3 + z/3 = 3. You are right! There are many solutions to this problem!

One adds more constraints and hope for unique solution. Give me a solution which has minimum length (minimize L₂-norm)? Answer is 3,3,3 (Generalized Inverse). Give me a solution which is most sparse (least number of non-zero entries)? Answer is 9,0,0 (or 0,0,9).

Many signals are sparse (in some domain). Images are sparse in Fourier/Wavelets domain; neural activity is sparse in time-domain, activity of programming is sparse in spatial-domain. Can we solve under-determined Ax = b for sparse x. Or formally,

min_xL₀(x), st Ax = b

where L₀ is the 0-norm i.e. number of non-zero entries.

This problem is hard to solve. L₀ is not a convex-function. Probably one can use convex relaxation trick and replace the non-convex function with convex one.

min_xL₁(x), st Ax = b

It is now well-understood that under certain assumptions on A, we can recover  by solving this tractable convex-optimization problem rather than the original NP hard problem. Greedy strategy based solvers works even faster though they are not proven to be as accurate as this one.

1 Somewhat formal statement

Skip it if formatting is not good.

If x_N (size N) is S − sparse then an under-determined system A_kNx_N = b_k -- where b is vector of k observations -- can be solved exactly given that A (measurement matrix) has following properties:

• ‖y‖² ∼ =‖x‖² where y = Ax for any sparse signal x i.e application of measurements matrix does not change the length of sparse signal x 'very much'. See Restricted Isometric Property in sec. 3.1 .

Lets there is a sparse signal x of length N. Its sparse. We take a dot-product of x with a random vector A_i of size N i.e. b_i = A_i * x. We make k such measurements. We can put all k A_i into a matrix A and represent the process by a linear system:

$\begin{bmatrix} A_1 \\ A_2 \\ \vdots \\ A_k \end{bmatrix} * x = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_k \end{bmatrix} \qquad(1)$

Note that each A_i has the dimension 1 × N. We can rewrite eq. 1 as following:

Ax = b  (2)
where dimension of A are k × N and k < <N. This is an under-determined system with the condition that x is sparse. Can we recover x (of size N) from b (of size k < < N)?

Figure fig. 1 shows a real system. Compressed sensing says that x can be recovered by solving the following linear program.

min_x‖x‖₁ given Ax = b  (3)

2 A demonstration using Python

Data for fig. 2 is generated by script ./compressed_sensing.py and for fig. 3, data is generated by script ./magic_reconstruction.py. Both of these scripts also generate similar figures.

Dependencies In addition to scipy, numpy, and matplotlib, we also need pyCSalgo. It is available via pip: pip install pyCSalgo --user.

Code is available at github

2.1 Algorithm

1. input : sparse signal x of size N
2. Generate a random matrix (measurement matrix) A of size k × N.
3. Make k measurements of x using A. That is b = Ax. Note that each measurement of x (i.e. entry of b) is some linear combination of values of x given by A_ix where A_i is ith row.
4. Now find x by solving min_x‖x‖₁given Ax = b.

For step 4, we are using function l1eq_pd from Python library pyCSalgo. This is re-implementation of l1magic routine written in Matlab by Candes. When k > >2S where S is the sparsity of x, we can recover x quite well. In practice, k ≥ 4S almost always gives exact result.

2.2 Demo 1: Sparse in time domain

Figure 1: x is the sparse signal of length 256. Using measurements matrix A, we made 64 random measurements of x namely b. The sparse solution of eq. 2 is the solution of min_x‖x‖₁given Ax = b.

Solution

Figure 2: Solution by CS which is almost equal to the original signal.

2.3 Sparse in Fourier domain

Figure 3: CS solution when signal is sparse in Fourier domain (DCT: Discrete Cosine Domain). Note that we took only 200 samples for a signal of size 2000. The recovery is pretty good. You can play both signal on phone and wouldn't notice a difference.

Appendix: Restricted isometric property (RIP)

Restricted isometric property (RIP) of a matrix A is defined as the following.

Given δ_s ∈ (0, 1), for any sub-matrix A_s of A, and for any sparse vector y, if following holds

(1 − δ_s)‖y‖² ≤ ‖A_sy‖² ≤ (1 + δ_s)‖y‖²
then matrix A is said to satify s−restricted isometric property with restricted isometric constant δ_s. Note that a matrix A with such a property does not change the length of singal x 'very much'. This enables us to sense two different sparse signal x₁ and x₂ such that Ax₁ and Ax₂ are almost likely to be different.

References

To learn about compressed sensing, I recommend following articles

• A very good motivational article www.ams.org/samplings/math-history/hap7-pixel.pdf . (last accessed Wed 30 Aug 2017 11:14:35 AM IST)
• A nice tutorial using Matlab CodeProject article.
• Terry Tao Blog with tag Compressed sensing.

There are many many other great articles and papers on this subject. There are some dedicated web-pages also.

- Dilawar

[Dilawar Singh]

Here is another demo : http://www.pyrunner.com/weblog/2016/05/26/compressed-sensing-python/

It has demo of images and uses cvxpy library. It is also much better written and exhaustive.

- Dilawar

[Kumar Ayush]

Thank you Dilawar. I hope you don't mind if I add your content to
wncc-iitb.org/wiki

Even being long, it is a good start for people wanting to learn about this

Regards
Ayush

--
--
The website for the club is http://wncc-iitb.org/
To post to this group, send email to wncc...@googlegroups.com
---
You received this message because you are subscribed to the Google Groups "Web and Coding Club IIT Bombay" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wncc_iitb+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Dilawar Singh]

On Thursday, September 7, 2017 at 10:33:32 AM UTC+5:30, cheekujodhpur wrote:

Thank you Dilawar. I hope you don't mind if I add your content to
wncc-iitb.org/wiki

Sure. Go ahead.