Practical Compressive Sensing with Toeplitz and Circulant Matrices

::News      ::Courses      ::People      ::Papers and Reports      ::Software

Compressive sensing encodes a signal into a relatively small number of incoherent linear measurements. In theory, the optimal incoherence is achieved by completely random measurement matrices. However, such matrices are difficult and/or costly to implement in hardware realizations.

Random Toeplitz and circulant matrices can be easily or even naturally realized in various applications. We introduce fast algorithms for reconstructing signals from incomplete Toeplitz and circulant measurements. Our computational results show that Toeplitz and circulant matrices are not only as effective as random matrices for signal encoding, but also permit much faster signal decoding.

Toeplitz and Circulant Matrices

Toeplitz and circulant matrices [wikipedia] have the forms, respectively, where every left-to-right descending diagonal is constant, i.e., Ti,j = Ti+1,j+1. If T satisfies the additional property that ti = tn+i, it is also a circulant matrix C.

Random Circulant vs I.I.D. Random Measurements

An m-by-n general matrix has mn degrees of freedom, but a circulant matrix of the same size has at most n degrees of freedom. Hence, a random circulant matrix is generated from much fewer independent random numbers or is "much less random" than an i.i.d. random matrix of the same size. This fact seemingly suggests that a random circulant matrix would yield less incoherent projections, and consequently worse CS recovery. However, numerical results below will hopefully convince us with a different story: circulant matrices can be equally effective as i.i.d. random matrices.

To obtain a broad picture, we test a variety of circulant matrix types and signal types for CS recovery. Each sensing matrix is generated by \$A=PC\$ from an n-by-n square circulant matrix \$C\$ and an m-by-n selection matrix \$P\$. We test seven different types of signals listed in the table above. A type-X3 signal can be written as \$x=\Psi u\$, where \$\Psi\$ is a discrete cosine transform (DCT) and \$u\$ is of type X1. Signals of type X6 often arise in signal and image processing. In all tests, we set the number of measurements to m=128 and signal dimension to n=512. For type X5, we fix sparsity k=8 and vary noise standard deviations sigma=1e-4*2i, for i=0,1,...,8. For each test, 1000 independent tests are run. Types X1 through X6 are recovered by L1 minimization with constraints Ax=b and type X7 by total variation minimization with the same constraints. Since the involved basis pursuit problems are small-scale linear programs, we use the highly reliable commercial package Gurobi version 2 with default settings. The comparison results given in Figures 1 show that on all except type-X3 signals, all the tested circulant matrices perform similarly and are no worse than i.i.d. random matrices (types I1 and I2) in terms of recovery frequencies and recovery errors. On the type-X3 signals, only matrices of types D1/P3 and D3/P3 can be competitive with the i.i.d. matrices; others give worse recovery results to different degrees. Although not all results are shown in the figures due to space limitation, performance differences have only been observed on the type-X3 signals.

To effectively encode type-X3 signals, it is important (see Romberg's 09 paper) to have uniform magnitude for \$di\$'s and random sampling (type P3). Besides, scrambling the orders of signal coordinates also works. All tested sensing matrices \$A\$ have the same performance on signals of type X4, which are obtained by randomly permuting the entries of type-X3 signals. However, the physical implementation of random permutations may become an added burden.

Fast algorithms and 2D results

Using Toeplitz and circulant sensing matrices allows significantly faster CS reconstruction compared to using i.i.d. random matrices. We introduce fast algorithms [download] for solving a variety of CS reconstruction models with Toeplitz and circulant sensing matrices, including equality fidelity, L1 and L2 square penalized fidelity, as well as one of or both L1 and total variation regularizations. They are portable to embedded and parallelized GPU computation. Furthermore, the algorithms can be extened to reconstruct three or higher-dimensional signals, either real or complex, and having single or multiple channels.

Partial Fourier sensing is a popular choice for sensing 2D images such as MR images. We compare images recovered from the same number of partial Fourier samples and partial circulant samples in Figure 2 below. At the same sample ratio of 11.21%, random circulant samples give better recovered images (Fig.2 (c) and (f)) than those from Fourier samples (Fig.2 (a), (b), (d), and (e)). Images (c) and (f) do not suffer from any aliasing artifacts of (a) and (d) or the loss of contrast of (b) and (e). Besides total variation reconstruction, we plan to adapt the proposed algorithms to enhanced reconstruction algorithms such as Bregman iterative regularization [link], nonconvex Lp (p<1) reconstrunction [pdf], and iterative support detection [pdf], and we expect to obtain improved results.

Here is a color image of the Great Wall of China reconstructed from 30% partial circulant measurements. Detecting dominant orientations in images without reconstruction It is rather difficult to learn geometric information of a signal directly from its CS samples without first recovering the signal, in general, because random CS samples are geometrically uncorrelated. We demonstrate that, however, it is possible to estimate edge orientations of images from partial circulant samples obtained from translating a rigid random mask. Figure 4 depicts two 192-by-192 images with different dominant edge orientations: the dominant orientation of the left image is horizontal, while that of the right image is vertical. Recall that 2D circulant samples are convolutions between the image and a moving mask (under circular boundary conditions). The measures of a random mask moving horizontally over the left image vary much less than those of the same mask moving vertically. We used this principle to decide whether an image being sensed was the one on the left or on the right of figure 4. For demonstrative purposes, we moved a mask 2s steps of one pixel each, to take totally m=2s+1 measures. The 2s steps form an L shape, out of which s are on the horizontal arm of L and another s are on its vertical arm. The entries of the mask took values from i.i.d. uniformly distributed random numbers in (0,1). We compared the TV value of all vertical-step measures to that of the horizontal-step measures to determine a dominant orientation according to the smaller TV value. The success rates were 97% for m=9, 99.7% for m=19, and 100% for m=49 out of 10000 independent trials for each m value.

Authors

Wotao Yin [Rice web], Simon Morgan (Los Alamos National Lab), Junfeng Yang [Nanjing web], Yin Zhang [Rice web]

Paper, Code, and References

• VCIP'10 Paper: [pdf]
W. Yin, S. Morgan, J. Yang, and Y. Zhang,
Practical Compressive Sensing with Toeplitz and Circulant Matrices,
In proceedings of Visual Communications and Image Processing (VCIP), 2010.
• MATLAB code: please email wotao dot yin @rice.edu for a free copy.
• References are included the VCIP'10 Paper above.

Ackowledgements

We thank Rick Chartrand (Los Alamos), Bingsheng He (Nanjing U.), Kevin Kelly (Rice) for valuable discussions.

The work of W. Yin was supported in part by NSF CAREER Award DMS-07-48839, ONR Grant N00014-08-1-1101, the U. S. Army Research Laboratory and the U. S. Army Research Office grant W911NF-09-1-0383, and an Alfred P. Sloan Research Fellowship. The work of S. Morgan was supported by the DOE. The work of J. Yang was supported in part by the Natural Science Foundation of China grant NSFC-10971095 and the Natural Science Foundation of Jiangsu Province BK2008255. The work of Y. Zhang was supported in part by NSF DMS-0811188 and ONR grant N00014-08-1-1101.