Compressed Sensing: A Python demo

I am having trouble with HTML format of blog post. Here is draft PDF PDF.

I recommend . This is much better written.

Can you solve $latex Ax=b$ for $latex 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 underdetermined 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!

If one adds more containts then one can hope to get a unique solution. Give me a solution which has minimum length (minimize (L_2)-norm)? Answer is 3,3,3 (hopefully). 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 acitivity is sparse in time-domain, activity of programming is sparse in spatial-domain. Can we solve under-determined (Ax=b) for sparse (x). Formally,

\min_x L_0(x), \quad \text{st} \quad Ax = b where (L_0) is the 0-norm i.e. number of non-zero entries.

This problem is hard to solve. (L_0) is not a convex function. In this case, one uses convex relaxation trick and replace the non-convex function with convex one.

\min_x L_1(x), \quad \text{st} \quad Ax = b

It is now well-understood that under certain assumptions on (A), we can recover (x) by solving this tractable convex-optimization problem rather than the original NP hard problem. Greedy strategy based solvers also exists which works even faster (though, according to Terry Tao blog, they do not guarantee the same performance as interior-point method based).

Somewhat formal statement

If (x_N) (size (N)) is (S-sparse) then an under-determined system (A_{kN}x_{N} =b_k) — where (b) is vector of (k) observations — can be solved exactly given that (A) (mesurement matrix) has following properties:

  • \left\Vert y\right\Vert^2 \sim = \left\Vert x \right\Vert ^2 where \(y = A x\) 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 \times N). We can rewrite eq. 1 as following:

A x = b\qquad(2) where dimension of (A) are (k \times 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 \left\Vert x \right\Vert_1 \; \text{given} \; Ax = b\qquad(3)

A demonstration using Python

Data for fig. 2 is generated by script ./ and for fig. 3, data is generated by script ./ 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


  1. input : sparse signal x of size N
  2. Generate a random matrix (measurement matrix) \(A\) of size \(k\times N\).
  3. Make \(k\) measurements of \(x\) using \(A\). That is \(b = A x\). Note that each measurement of \(x\) (i.e. entry of \(b\)) is some linear combination of values of \(x\) given by \(A_i x\) where \(A_i\) is ith row.
  4. Now find \(x\) by solving \(\min_x \left\Vert x \right\Vert_1 \text{given}\; Ax=b\).

For step 4, we are using function l1eq_pd from Python library pyCSalgo. This is reimplementation 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 \ge 4S) almost always gives exact result.

There are other algorithms available which works faster; they are based on greedy strategy. We are not discussing them here.

Sparse in time domain


Figure 1: x is the sparse signal of length 256. Using measurements matrix A, we made 124 random measurements of \(x\) namely b. The sparse solution of eq. 2 is the solution of $\min_x \left\Vert x \right\Vert_1 \text{given}\; Ax =b$.

Sparse in Fourier domain


Figure 3: CS solution when signal is sparse in Forier domain. Note that we took 200 samples for a singal of size 2000. The recovery is pretty good.

With error in measurements

This example shows that some error is added to measurements. The reconstructed signal is no longer ‘exact’ enough. But still, given that x was binary signal, we can recover x by thresholding function.


Mathematical definitions

Restricted isometric property (RIP)

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

Given (\delta_s \in (0,1)), for any sub-matrix (A_s) of (A), and for any sparse vector (y), if following holds

(1-\delta_s) \left\Vert y \right\Vert^2 \leq \left\Vert A_s y \right\Vert ^2 \leq (1+\delta_s) \left\Vert y \right\Vert ^2

then matrix (A) is said to satisfy (s-)restricted isometric property with restricted isometric constant (\delta_s). Note that a matrix A with such a property does not change the length of signal (x) ‘very much’. This enables us to sense two different sparse signal (x_1) and (x_2) such that (A x_1) and (A x_2) are almost likely to be different.


To learn about compressed sensing, I recommend following articles

There are many many other great articles and papers on this subject. There are some dedicated webpages also.


Performance of “sorting dictionary by values” in python2, python3 and pypy

The script is hosted here . It is based on the work of

My script has been changed to accommodate python3 (iteritems is gone and replaced by items — not sure whether it is a fair replacement). For method names and how they are implemented, please refer to script or the blog post.

Following chart shows the comparison. PyPy does not boost up the performance for simple reason that dictionary sorted is not large enough. I’ve put it here just for making a point and PyPy can slow thing down on small size computation.

The fastest method is sbv6 which is based on PEP-0265 is the fastest. Python3 always performing better than python2.




Writing Maxima expression to text file in TeX format (for LaTeX)

You want to write an Maxima expression to a file which can be read by other application e.g. LaTeX.

Lets say the expression is sys which contains variable RM. You first want to replace RM by R_m .  Be sure to load mactex-utilities if you have matrix. Without loading this module, the tex command generates TeX output, not LaTeX.

load( "mactex-utilities" )$ 
sys : RM * a / b * log( 10 )$
texput( RM, "R_m")$
sysTex : tex( sys, false)$
with_stdout( "outout.txt", display( sysTex ) )$

Other methods such as stringout, save and write put extra non-TeX characters in file.

I get the following in file outout.txt after executing the above.

{{\log 10\,R_m\,a}\over{b}}

Image stabilization using OpenCV

This application deals with video of neural recordings. In such recordings, feature sizes are small. On top of it, recordings are quite noisy. Animal head movements introduces sharp shakes. Out of the box video stabilizer may not work very well on such recordings. Though there are quite a lot of plugins for ImageJ to do such a work, I haven’t compared their performance with this application. This application is hosted here and a demo video is available on youtube here .

The summary of basic principle is following:

0. Collect all frames in a list/vector.

1. Use bilateral filter to smooth out each frame. Bilateral filter smoothens image without distorting the edges (well to a certain extent).

2.  Calculate optical flow between previous frame and current frame. This is a proxy for movement. Construct a transformation and store them in a vector. OpenCV function `goodFeatureToTrack` does almost all the work for us.

  1. Take average of these transformations and apply it on each frame of original recording; that’s correct motion.

A csv reader based on Haskell-cassava library : performance

I implemented my own csv reader using cassava library. The reader from missingh library was taking too long (~ 17 seconds) for a file with 43200 lines. I compared the result with python-numpy and python-pandas csv reader. Below is rough comparison.

cassava (ignore #) 3.3 sec
cassava (no support for ignoring #) 2.7 sec
numpy loadtxt > 10 sec
pandas read_csv 1.5 sec

As obvious, pandas does really well at reading csv file. I was hoping that my csv reader would do better but it didn’t. But it still beats the parsec based reader hands down.

The code is here

Thresholding numpy array, well sort of

Here is a test case

>>> import numpy as np
>>> a = np.array( [ 0.0, 1, 2, 0.2, 0.0, 0.0, 2, 3] )

I want to turn all non-zero elements of this array to 1. I can do it using np.where and numpy indexing.

>>>  a[ np.where( a != 0 ) ] = 1
>>> a
array([ 0.,  1.,  1.,  1.,  0.,  0.,  1.,  1.])

note: np.where returns the indices where the condition is true. e.g. if you want to change all 0 to -1.

>>> a[ np.where[ a == 0] ] = -1.0

That’s it. Checkout np.clip as well.