Lasso Optimisation

This example demonstrates the use of class bpdn.BPDNProjL1 to solve the least absolute shrinkage and selection operator (lasso) problem [28]

\[\mathrm{argmin}_\mathbf{x} \; (1/2) \| D \mathbf{x} - \mathbf{s} \|_2^2 \; \text{such that} \; \| \mathbf{x} \|_1 \leq \gamma\]

where \(D\) is the dictionary, \(\mathbf{x}\) is the sparse representation, and \(\mathbf{s}\) is the signal to be represented. In this example the lasso 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.admm import bpdn
from sporco import plot

Configure problem size, sparsity, and noise level.

N = 512      # Signal size
M = 4*N      # Dictionary size
L = 32       # Number of non-zero coefficients in generator
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
D = np.random.randn(N, M)
x0 = np.zeros((M, 1))
si = np.random.permutation(list(range(0, M-1)))
x0[si[0:L]] = np.random.randn(L, 1)

# Construct reference and noisy signal
s0 =
s = s0 + sigma*np.random.randn(N,1)

Set bpdn.BPDNProjL1 solver class options. The value of \(\gamma\) has been manually chosen for good performance.

gamma = 5.0
opt = bpdn.BPDNProjL1.Options({'Verbose': True, 'MaxMainIter': 500,
                    'RelStopTol': 1e-6, 'AutoRho': {'RsdlTarget': 1.0}})

Initialise and run BPDNProjL1 object

b = bpdn.BPDNProjL1(D, s, gamma, opt)
x = b.solve()

print("BPDNProjL1 solve time: %.2fs" % b.timer.elapsed('solve'))
Itn   Fnc       Cnstr     r         s         ρ
   0  4.85e+03  0.00e+00  8.84e-01  3.18e-01  1.00e+00
   1  4.55e+03  0.00e+00  6.20e-01  1.41e-01  1.00e+00
   2  4.74e+03  0.00e+00  6.34e-01  5.22e-02  1.00e+00
   3  4.66e+03  0.00e+00  5.58e-01  1.52e-02  1.00e+00
   4  4.53e+03  0.00e+00  6.29e-01  3.80e-02  1.00e+00
   5  4.60e+03  0.00e+00  6.22e-01  2.88e-02  1.00e+00
   6  4.61e+03  0.00e+00  6.75e-01  1.02e-02  1.00e+00
   7  4.55e+03  0.00e+00  6.30e-01  1.75e-02  1.00e+00
   8  4.59e+03  2.72e-16  6.25e-01  1.74e-02  1.00e+00
   9  4.61e+03  1.03e-15  5.79e-01  9.59e-03  1.00e+00
  10  4.54e+03  0.00e+00  5.85e-01  3.14e-02  7.77e+00
  11  4.51e+03  6.28e-16  6.06e-01  4.06e-02  7.77e+00
  12  4.55e+03  6.28e-16  6.25e-01  2.00e-02  7.77e+00
  13  4.53e+03  6.28e-16  6.26e-01  7.87e-03  7.77e+00
  14  4.52e+03  2.72e-16  6.12e-01  1.67e-02  7.77e+00
  15  4.55e+03  0.00e+00  5.88e-01  1.25e-02  7.77e+00
  16  4.54e+03  0.00e+00  5.74e-01  5.31e-03  7.77e+00
  17  4.51e+03  0.00e+00  5.78e-01  8.35e-03  7.77e+00
  18  4.51e+03  0.00e+00  5.91e-01  8.02e-03  7.77e+00
  19  4.52e+03  2.24e-15  5.99e-01  4.72e-03  7.77e+00
  20  4.51e+03  6.28e-16  5.81e-01  1.84e-02  8.76e+01
  21  4.51e+03  0.00e+00  5.43e-01  2.32e-02  8.76e+01
  22  4.51e+03  2.72e-16  5.04e-01  1.15e-02  8.76e+01
  23  4.51e+03  0.00e+00  4.75e-01  4.53e-03  8.76e+01
  24  4.50e+03  0.00e+00  4.56e-01  8.82e-03  8.76e+01
  25  4.50e+03  0.00e+00  4.41e-01  6.69e-03  8.76e+01
  26  4.50e+03  6.28e-16  4.23e-01  2.56e-03  8.76e+01
  27  4.50e+03  0.00e+00  4.00e-01  3.83e-03  8.76e+01
  28  4.50e+03  1.26e-15  3.75e-01  3.85e-03  8.76e+01
  29  4.50e+03  0.00e+00  3.52e-01  2.24e-03  8.76e+01
  30  4.50e+03  3.51e-16  2.54e-01  8.91e-03  1.10e+03
  31  4.49e+03  0.00e+00  1.34e-01  1.38e-02  1.10e+03
  32  4.49e+03  0.00e+00  7.39e-02  1.13e-02  1.10e+03
  33  4.49e+03  6.28e-16  4.18e-02  7.15e-03  1.10e+03
  34  4.49e+03  3.68e-16  2.39e-02  3.56e-03  1.10e+03
  35  4.49e+03  0.00e+00  1.38e-02  1.38e-03  1.10e+03
  36  4.49e+03  0.00e+00  8.10e-03  5.96e-04  1.10e+03
  37  4.49e+03  0.00e+00  4.81e-03  4.61e-04  1.10e+03
  38  4.49e+03  0.00e+00  2.90e-03  3.04e-04  1.10e+03
  39  4.49e+03  3.68e-16  1.77e-03  1.46e-04  1.10e+03
  40  4.49e+03  2.94e-16  7.08e-04  8.66e-05  3.82e+03
  41  4.49e+03  6.28e-16  1.04e-04  1.23e-04  3.82e+03
  42  4.49e+03  0.00e+00  3.05e-05  9.93e-05  3.82e+03
  43  4.49e+03  0.00e+00  2.14e-05  7.83e-05  3.82e+03
  44  4.49e+03  0.00e+00  1.63e-05  6.01e-05  3.82e+03
  45  4.49e+03  0.00e+00  1.26e-05  4.63e-05  3.82e+03
  46  4.49e+03  0.00e+00  9.65e-06  3.55e-05  3.82e+03
  47  4.49e+03  0.00e+00  7.43e-06  2.73e-05  3.82e+03
  48  4.49e+03  0.00e+00  5.71e-06  2.10e-05  3.82e+03
  49  4.49e+03  0.00e+00  4.40e-06  1.62e-05  3.82e+03
  50  4.49e+03  0.00e+00  6.44e-06  1.22e-05  1.99e+03
  51  4.49e+03  0.00e+00  4.88e-06  7.99e-06  1.99e+03
  52  4.49e+03  0.00e+00  2.99e-06  4.48e-06  1.99e+03
  53  4.49e+03  0.00e+00  1.69e-06  2.34e-06  1.99e+03
  54  4.49e+03  0.00e+00  8.96e-07  1.16e-06  1.99e+03
  55  4.49e+03  0.00e+00  4.56e-07  5.61e-07  1.99e+03
BPDNProjL1 solve time: 0.12s

Plot comparison of reference and recovered representations.

plot.plot(np.hstack((x0, x)), title='Sparse representation',
          lgnd=['Reference', 'Reconstructed'])

Plot functional value, residuals, and rho

its = b.getitstat()
fig = plot.figure(figsize=(20, 5))
plot.subplot(1, 3, 1)
plot.plot(its.ObjFun, xlbl='Iterations', ylbl='Functional', fig=fig)
plot.subplot(1, 3, 2)
plot.plot(np.vstack((its.PrimalRsdl, its.DualRsdl)).T,
          ptyp='semilogy', xlbl='Iterations', ylbl='Residual',
          lgnd=['Primal', 'Dual'], fig=fig)
plot.subplot(1, 3, 3)
plot.plot(its.Rho, xlbl='Iterations', ylbl='Penalty Parameter', fig=fig)