CVXMOD example – Spike signal reconstruction

Given a measurement matrix

and a data signal

, this code minimizes |y - Ax|_2 + gamma |x|_p

for

and

and produces the below figure. Note how l_1

regularized reconstruction preserves spikes in the signal, while l_2

regularization would heavily penalize such large deviations.

"""Signal reconstruction example with various norms.

Written by Kwangmoo Koh, 2007-03-09.
Adapted for CVXMOD by Timothy Hunter, 2008-05-10.

This file shows how to reconstruct a spike signal using l1 or l2
regularization."""

from cvxmod import *
import pylab

def l1reg(A, y, gamma):
    """Given data y and measurement matrix A, returns a reconstructed signal
    using an l1 heuristic."""
    n = cols(A)
    x1 = optvar('x1', n)
    p = problem(minimize(norm2(y - A*x1) + gamma*norm1(x1)))
    p.solve()
    return value(x1)

def l2reg(A, y, gamma):
    """Given data y and measurement matrix A, returns a reconstructed signal
    using an l2 heuristic."""
    n = cols(A)
    x1 = optvar('x1', n)
    p = problem(minimize(norm2(y - A*x1) + gamma*norm2(x1)))
    p.solve()
    return value(x1)

# Analyze a particular signal.
m = 50 # number of observations.
n = 200 # signal length.
s = 1.5e-2 # spike density.

# Generate a spike signal.
randseed(5)
x0 = zeros(n, 1)
for i in range(n):
    if rand() < s:
        x0[i] = 1

# Generate a random measurement matrix.
A = randn(m, n)

# Generate noisy measurements.
sigma = 0.1
v = sigma*randn(m, 1)
y = A*x0 + v

# Find signal reconstructions.
gamma = 0.01
x1 = l1reg(A, y, gamma)
x2 = l2reg(A, y, gamma)

if True: # should we display figures?
    axis = [0, n, -1.2, 1.2]

    pylab.subplot(311)
    pylab.stem(range(n), matrix(x0), 'k-', 'k.', 'k')
    pylab.axis(axis)
    pylab.ylabel('original')

    pylab.subplot(312)
    pylab.stem(range(n), matrix(x1), 'k-', 'k.', 'k')
    pylab.axis(axis)
    pylab.ylabel('l1 regularized')

    pylab.subplot(313)
    pylab.stem(range(n), x2, 'k-', 'k.', 'k')
    pylab.axis(axis)
    pylab.ylabel('l2 regularized')
    pylab.subplots_adjust(0.08, 0.04, 0.98, 0.98)
    pylab.savefig('www/examples/figures/spike.png', dpi=65)
show()

CVXMOD example – Spike signal reconstruction

Reconstructed signals

Python and CVXMOD code (spike.py)