# Basis Pursuit DeNoising with Joint Sparsity¶

This example demonstrates the use of class bpdn.BPDNJoint to solve the Basis Pursuit DeNoising (BPDN) problem with joint sparsity via an $$ℓ_{2,1}$$ norm term

$\mathrm{argmin}_\mathbf{x} \; (1/2) \| D X - S \|_2^2 + \lambda \| X \|_1 + \mu \| X \|_{2,1}$

where $$D$$ is the dictionary, $$X$$ is the sparse representation, and $$S$$ is the signal to be represented. In this example the BPDN problem is used to estimate the reference sparse representation that generated a signal from a noisy version of the signal.

from __future__ import print_function
from builtins import input
from builtins import range

import numpy as np

from sporco import util
from sporco import plot
plot.config_notebook_plotting()


Configure problem size, sparsity, and noise level.

N = 32         # Signal size
M = 4*N        # Dictionary size
L = 12         # Number of non-zero coefficients in generator
K = 16         # Number of signals
sigma = 0.5    # Noise level


Construct random dictionary, reference random sparse representation, and test signal consisting of the synthesis of the reference sparse representation with additive Gaussian noise.

# Construct random dictionary and random sparse coefficients
np.random.seed(12345)
D = np.random.randn(N, M)
x0 = np.zeros((M, K))
si = np.random.permutation(list(range(0, M-1)))
x0[si[0:L],:] = np.random.randn(L, K)
# Construct reference and noisy signal
s0 = D.dot(x0)
s = s0 + sigma*np.random.randn(N,K)


Set BPDNJoint solver class options.

opt = bpdn.BPDNJoint.Options({'Verbose': False, 'MaxMainIter': 500,
'RelStopTol': 1e-3, 'rho': 10.0,
'AutoRho': {'RsdlTarget': 1.0}})


Select regularization parameters $$\lambda, \mu$$ by evaluating the error in recovering the sparse representation over a logarithmicaly spaced grid. (The reference representation is assumed to be known, which is not realistic in a real application.) A function is defined that evalues the BPDN recovery error for a specified $$\lambda, \mu$$, and this function is evaluated in parallel by sporco.util.grid_search.

# Function computing reconstruction error for (lmbda, mu) pair
def evalerr(prm):
lmbda = prm[0]
mu = prm[1]
b = bpdn.BPDNJoint(D, s, lmbda, mu, opt)
x = b.solve()
return np.sum(np.abs(x-x0))

# Parallel evalution of error function on lmbda,mu grid
lrng = np.logspace(-4, 0.5, 10)
mrng = np.logspace(0.5, 1.6, 10)
sprm, sfvl, fvmx, sidx = util.grid_search(evalerr, (lrng, mrng))
lmbda = sprm[0]
mu = sprm[1]

print('Minimum ℓ1 error: %5.2f at (𝜆,μ) = (%.2e, %.2e)' % (sfvl, lmbda, mu))

Minimum ℓ1 error: 40.36 at (𝜆,μ) = (1.00e-04, 1.71e+01)


Once the best $$\lambda, \mu$$ have been determined, run bpdn.BPDNJoint.solve with verbose display of ADMM iteration statistics.

# Initialise and run BPDNJoint object for best lmbda and mu
opt['Verbose'] = True
b = bpdn.BPDNJoint(D, s, lmbda, mu, opt)
x = b.solve()

print("BPDNJoint solve time: %.2fs" % b.timer.elapsed('solve'))

Itn   Fnc       DFid      Regℓ1     Regℓ2,1   r         s         ρ
--------------------------------------------------------------------------
0  2.50e+03  2.44e+03  1.08e+01  3.18e+00  9.14e-01  1.43e-01  1.00e+01
1  1.01e+03  5.79e+02  8.48e+01  2.55e+01  5.50e-01  3.77e-01  1.00e+01
2  1.20e+03  4.44e+02  1.46e+02  4.39e+01  2.81e-01  3.51e-01  1.00e+01
3  8.65e+02  1.72e+02  1.35e+02  4.05e+01  1.99e-01  2.26e-01  1.00e+01
4  7.83e+02  1.76e+02  1.19e+02  3.55e+01  1.51e-01  1.96e-01  1.00e+01
5  7.74e+02  2.04e+02  1.13e+02  3.33e+01  1.27e-01  1.00e-01  1.00e+01
6  7.39e+02  1.54e+02  1.16e+02  3.42e+01  9.41e-02  7.47e-02  1.00e+01
7  7.30e+02  1.08e+02  1.23e+02  3.63e+01  6.86e-02  7.00e-02  1.00e+01
8  7.27e+02  9.87e+01  1.24e+02  3.67e+01  5.50e-02  4.15e-02  1.00e+01
9  7.23e+02  1.12e+02  1.20e+02  3.57e+01  4.29e-02  3.56e-02  1.00e+01
10  7.22e+02  1.21e+02  1.18e+02  3.51e+01  3.34e-02  2.73e-02  1.10e+01
11  7.20e+02  1.09e+02  1.20e+02  3.57e+01  2.62e-02  2.32e-02  1.10e+01
12  7.19e+02  9.89e+01  1.22e+02  3.63e+01  2.12e-02  1.72e-02  1.10e+01
13  7.19e+02  1.01e+02  1.21e+02  3.61e+01  1.65e-02  1.05e-02  1.10e+01
14  7.19e+02  1.05e+02  1.21e+02  3.59e+01  1.31e-02  9.99e-03  1.10e+01
15  7.19e+02  1.05e+02  1.21e+02  3.59e+01  1.07e-02  7.80e-03  1.10e+01
16  7.19e+02  1.01e+02  1.21e+02  3.61e+01  8.85e-03  5.53e-03  1.10e+01
17  7.18e+02  1.01e+02  1.21e+02  3.61e+01  7.16e-03  3.25e-03  1.10e+01
18  7.18e+02  1.02e+02  1.21e+02  3.60e+01  5.79e-03  3.35e-03  1.10e+01
19  7.18e+02  1.02e+02  1.21e+02  3.60e+01  4.95e-03  2.53e-03  1.10e+01
20  7.18e+02  1.01e+02  1.21e+02  3.61e+01  3.93e-03  1.73e-03  1.54e+01
21  7.18e+02  1.01e+02  1.21e+02  3.61e+01  3.18e-03  1.78e-03  1.54e+01
22  7.18e+02  1.01e+02  1.21e+02  3.61e+01  2.63e-03  1.60e-03  1.54e+01
23  7.18e+02  1.01e+02  1.21e+02  3.61e+01  2.19e-03  1.07e-03  1.54e+01
24  7.18e+02  1.01e+02  1.21e+02  3.61e+01  1.81e-03  8.09e-04  1.54e+01
25  7.18e+02  1.01e+02  1.21e+02  3.61e+01  1.49e-03  7.34e-04  1.54e+01
26  7.18e+02  1.01e+02  1.21e+02  3.61e+01  1.24e-03  5.96e-04  1.54e+01
27  7.18e+02  1.01e+02  1.21e+02  3.61e+01  1.04e-03  4.27e-04  1.54e+01
28  7.18e+02  1.01e+02  1.21e+02  3.61e+01  8.69e-04  3.11e-04  1.54e+01
--------------------------------------------------------------------------
BPDNJoint solve time: 0.04s


Plot comparison of reference and recovered representations.

fig = plot.figure(figsize=(6, 8))
plot.subplot(1, 2, 1)
plot.imview(x0, cmap=plot.cm.Blues, title='Reference', fig=fig)
plot.subplot(1, 2, 2)
plot.imview(x, cmap=plot.cm.Blues, title='Reconstruction', fig=fig)
fig.show()


Plot lmbda,mu error surface, functional value, residuals, and rho

its = b.getitstat()
fig = plot.figure(figsize=(15, 10))
ax = fig.add_subplot(2, 2, 1, projection='3d')
ax.locator_params(nbins=5, axis='y')
plot.surf(fvmx, x=np.log10(mrng), y=np.log10(lrng), xlbl='log($\mu$)',
ylbl='log($\lambda$)', zlbl='Error', fig=fig, ax=ax)
plot.subplot(2, 2, 2)
plot.plot(its.ObjFun, xlbl='Iterations', ylbl='Functional', fig=fig)
plot.subplot(2, 2, 3)
plot.plot(np.vstack((its.PrimalRsdl, its.DualRsdl)).T,
ptyp='semilogy', xlbl='Iterations', ylbl='Residual',
lgnd=['Primal', 'Dual'], fig=fig)
plot.subplot(2, 2, 4)
plot.plot(its.Rho, xlbl='Iterations', ylbl='Penalty Parameter', fig=fig)
fig.show()