> ICTEAM > ELEN > ISPGroup > Laurent Jacques' Homepage > QIHTDemo

Abstract

"Quantized Iterative Hard Thresholding algorithm demo" associated to the paper:

"Quantized Iterative Hard Thresholding: Bridging 1-bit and High-Resolution Quantized Compressed Sensing"
Laurent Jacques, Kvin Degraux, Christophe De Vleeschouwer,
arXiv:1305.1786

Updates:

  • June 11, 2013: Initial release of this small demo.

Download the demo

You can donwload the single Matlab mfile in zip format here.


Published Demo

Contents

% DEMONSTRATION file for the paper:
%
% Quantized Iterative Hard Thresholding: Bridging 1-bit and High-Resolution
% Quantized Compressed Sensing
%
% Laurent Jacques, Kévin Degraux and Christophe De Vleeschouwer
% June 2013
% -------------------------------------------------------------------------
%
% This simple demo file provides a rapid overview of the Quantized
% Iterative Hard Thresholding (QIHT) [1] capability in the Quantized
% Compressed Sensing framework. This demo is realized by a comparison with
% a standard implementation of Iterative Hard Thresholding [2].
%
% Given a B-bits scalar quantizer Q() applied pointwise on
% vectors, let us take the common compressed sensing setup where a K-sparse
% signal x in R^N is sensed through the model
%
%       y = A*x in R^M,
%
% for some sensing matrix (e.g., random Gaussian) of size MxN.
%
% QIHT aims at estimating x from
%
%          q = Q(y) = Q(A*x),
%
% encoded over R = BM bits knowing K and A.
%
% As shown in [1], the simple fact of inserting Q in the Iterative Hard
% Thresholding algorithm, i.e., defining
%
%      x0 = 0;
%      xn+1 = H_K(xn + mu * A'*(q - Q(A*xn));
%
% for H_K the best K-term hard thresholding operator and mu>0, is
% theoretically valid. Indeed, at each iteration the inconsistency energy
%
%      E(q, A*u) = sum_i J((A*u)_i, q_i)
%
% is minimized with
%
%      J(a,b) = sum_{k=2}^{2^B-1} (w_{k+1}-w_{k}) B(a - tau_k, b - tau_k),
%
%      B(a,b) = |(sign(b)*a)_-| := .5*(|sign(b)*a| - (sign(b)*a)).
%
% where w_k and tau_k are the levels and thresholds of the scalar quantizer
% Q, respectively.
%
% J has two interesting limit cases. For B=1,
%
%         J(a,b) = 2w_0 |(sign(b)*a)_-|
%
% which is the energy minized by the Binary Iterative Hard Thresholding
% algorithm in [3] for estimating a sparse signal from its 1-bit CS
% measurements w0*sign(A*x).
%
% And for B>>1, we can show that J(a,b) -> .5 |a - b|^2, i.e., the
% quadratic energy minmized by IHT [2].
%
% In [1], the value mu = c/M, with c = (1 - sqrt(3*K/M)), is shown to give
% good convergence results for random Gaussian matrices.
%
% References:
%
% [1] Laurent Jacques, Kévin Degraux, Christophe De
% Vleeschouwer, "Quantized Iterative Hard Thresholding: Bridging 1-bit and
% High-Resolution Quantized Compressed Sensing"
% http://arxiv.org/abs/1305.1786
%
% [2] T. Blumensath and M.E. Davies, "Iterative hard thresholding for
% compressed sensing," App. Comp. Harm. Anal., 27(3):265-274, 2009.
%
% [3] L. Jacques, J.N. Laska, P.T. Boufounos, and R.G. Baraniuk, ?Robust
% 1-bit compressive sensing via binary stable embeddings of sparse
% vectors,? IEEE Trans. Inf. Th., 59(4):2082-2102, 2013.

First, let us set some dimensions ...

N = 1024;          % ambient signal space dimension
B = 2;             % number of bits per measurements
M = round(2048/B); % number of measurements
R = B*M;           % constant rate R = BM
K = 16;            % signal sparsity level (canonical basis)

nbins = 2^B;

% ... and some useful anonymous functions

% parenthesis selection
par = @(in, varargin) in(varargin{:});

% hard thresholding operator (valid for all distinct component)
% caveat: does not work on vectors with non-distinct heights
H = @(in,k) in.*(abs(in) >= par(sort(abs(in), 'descend'), k));

% Normalized SNR
unit = @(in) in(:)/norm(in(:));
NSNR = @(in1,in2) 20*log10(1/norm(unit(in1) - unit(in2)));

Generating a unit K-sparse signal.

Imposing unit norm allows one fair comparison with B=1 case (with total lost of amplitude)

x = zeros(N,1);
ind = randperm(N);
x(ind(1:K))=randn(K,1);
x = x/norm(x);

Random Gaussian Sensing Matrix

A = randn(M,N);

Lloyd-Max quantizer optimizer for normally distributed inputs

Reference: R. M. Gray and D. L. Neuhoff, "Quantization," IEEE Trans. Inf. Th., 44(6):2325-2383, 1998.

switch B
    case 1, % 2=2^1 bins
        thr = [-Inf, 0, Inf];
        lev = [-0.7979, 0.7979];
    case 2, % 4=2^2 bins
        thr = [-Inf   -0.9816    0.0000    0.9816       Inf];
        lev = [-1.5104   -0.4528    0.4528    1.5104];
    case 3, % 8=2^3 bins
        thr = [-Inf   -1.7479   -1.0500   -0.5005   -0.0000 ...
            0.5005    1.0500    1.7479       Inf];
        lev = [-2.1519   -1.3439   -0.7560   -0.2451    0.2451 ...
            0.7560    1.3439    2.1519];
    case 4, % 16=2^4 bins
        thr = [-Inf   -2.4008   -1.8435   -1.4371   -1.0993  ...
            -0.7995   -0.5224   -0.2582    0.0000    0.2582  ...
             0.5224    0.7995    1.0993    1.4371    1.8435  ...
             2.4008       Inf];
         lev = [-2.7326   -2.0690   -1.6180   -1.2562   -0.9423  ...
             -0.6568   -0.3880   -0.1284    0.1284    0.3880    0.6568 ...
             0.9423    1.2562    1.6180    2.0690    2.7326];
    otherwise
        error('This demo file works only for 1<=B<=4');
end


Qc = @(x) lev([x >= thr(1:end-1)] & [x < thr(2:end)]);
Q = @(x) arrayfun(Qc,x);

Quantized Compressed Sensing of x

y = A*x;
q = Q(y);

% showing the results
figure(1); clf; hold on;
title('Visually comparing unquantized (black) and quantized measurements (red)');
plot(y,'k-', 'linewidth', 2);
plot(q,'ro', 'linewidth', 2);
ylim([-3, 3]);
axis tight;
xlabel('measurement index')
for k=1:nbins,
plot(y*0+lev(k),'g--', 'linewidth', 2);
end

Reconstructing the signal with Iterative Hard Thresholding

xn = zeros(N,1);

c = 1 - sqrt(3*K/M);
mu = c/M;
maxiter = 1000;
stopcrit = 1e-7;

for n = 1:maxiter
    oxn = xn;

    % main IHT iteration
    a = xn + mu * A' * (q - A*xn);
    xn = H(a,K);

    if (norm(oxn-xn)/norm(xn) < stopcrit)
        break;
    end
end

fprintf('\nQCS dims: N:%i, K: %i, M:%i, B:%i, R=B*M: %i bits\n', N, K, M, B, B*M);
fprintf('IHT  -> Nb iter: %4.f, ||Ax-q||=%e, ||Q(Ax)-q||=%e, SNR:%3.2fdB\n', ...
    n, norm(A*xn - q), norm(Q(A*xn) - q), NSNR(x,xn) );

iht_x = xn;
QCS dims: N:1024, K: 16, M:1024, B:2, R=B*M: 2048 bits
IHT  -> Nb iter:   17, ||Ax-q||=1.011685e+01, ||Q(Ax)-q||=8.007934e+00, SNR:25.53dB

Reconstructing the signal with Quantized Iterative Hard Thresholding

xn = zeros(N,1);

c = 1 - sqrt(3*K/M);
mu = c/M;
maxiter = 1000;
stopcrit = 1e-7;

for n = 1:maxiter
    oxn = xn;


    % main QIHT iteration
    a = xn + mu * A' * (q - Q(A*xn));
    xn = H(a,K);


    if (norm(oxn-xn)/norm(xn) < stopcrit)
        break;
    end
end

fprintf('\nQCS dims: N:%i, K: %i, M:%i, B:%i, R=B*M: %i bits\n', N, K, M, B, B*M);
fprintf('QIHT -> Nb iter: %4.f, ||Ax-q||=%e, ||Q(Ax)-q||=%e, SNR:%3.2fdB\n', ...
    n, norm(A*xn - q), norm(Q(A*xn) - q), NSNR(x,xn) );

qiht_x = xn;
QCS dims: N:1024, K: 16, M:1024, B:2, R=B*M: 2048 bits
QIHT -> Nb iter:   23, ||Ax-q||=1.071389e+01, ||Q(Ax)-q||=0.000000e+00, SNR:45.10dB

Showing results

% signals
figure(2); clf;

subplot(3,1,1); hold on;
plot(unit(x),'linewidth',2);
title('Original signal x');
yl = 1.2*ylim;
ylim(yl);

subplot(3,1,2); hold on;
plot(unit(iht_x),'linewidth',2);
title(sprintf('IHT estimate, Normalized SNR: %2.2f dB   (x in red)', NSNR(x, iht_x)));
ylim(yl);
plot(find(x), x(x~=0), 'rx', 'linewidth', 2);

subplot(3,1,3); hold on;
plot(unit(qiht_x),'linewidth',2);
title(sprintf('QIHT estimate, Normalized SNR: %2.2f dB   (x in red)', NSNR(x, qiht_x)));
ylim(yl);
plot(find(x), x(x~=0), 'rx', 'linewidth', 2);

% differences
figure(3); clf;

subplot(2,1,1); hold on;
plot(unit(x) - unit(iht_x),'linewidth',2);
title(sprintf('IHT reconstruction error, max error: %e   (support of x in red)', norm(x - iht_x, Inf)));
plot(find(x), 0*x(x~=0), 'rx', 'linewidth', 2);


subplot(2,1,2); hold on;
plot(unit(x) - unit(qiht_x),'linewidth',2);
title(sprintf('QIHT reconstruction error, max error: %e   (support of x in red)', norm(x - qiht_x, Inf)));
plot(find(x), 0*x(x~=0), 'rx', 'linewidth', 2);