Basis Pursuit DeNoising (APGM)

This example demonstrates the solution of the the sparse coding problem

\[\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - D \mathbf{x} \|_2^2 + \lambda \| \mathbf{x} \|_1\;,\]

where \(D\) the dictionary, \(\mathbf{y}\) the signal to be represented, and \(\mathbf{x}\) is the sparse representation.

import numpy as np

import komplot as kplt

import scico.numpy as snp
from scico import functional, linop, loss
from scico.optimize.pgm import AcceleratedPGM
from scico.util import device_info
kplt.config_notebook_plotting()

Construct a random dictionary, a reference random sparse representation, and a test signal consisting of the synthesis of the reference sparse representation.

m = 512  # Signal size
n = 4 * m  # Dictionary size
s = 32  # Sparsity level (number of non-zeros)
σ = 0.5  # Noise level

np.random.seed(12345)
D = np.random.randn(m, n).astype(np.float32)
L0 = np.linalg.norm(D, 2) ** 2

x_gt = np.zeros(n, dtype=np.float32)  # true signal
idx = np.random.permutation(list(range(0, n - 1)))
x_gt[idx[0:s]] = np.random.randn(s)
y = D @ x_gt + σ * np.random.randn(m)  # synthetic signal

x_gt = snp.array(x_gt)  # convert to jax array
y = snp.array(y)  # convert to jax array

Set up the forward operator and AcceleratedPGM solver object.

maxiter = 100
λ = 2.98e1
A = linop.MatrixOperator(D)
f = loss.SquaredL2Loss(y=y, A=A)
g = λ * functional.L1Norm()
solver = AcceleratedPGM(
    f=f, g=g, L0=L0, x0=A.adj(y), maxiter=maxiter, itstat_options={"display": True, "period": 10}
)

Run the solver.

print(f"Solving on {device_info()}\n")
x = solver.solve()
hist = solver.itstat_object.history(transpose=True)

Solving on GPU (NVIDIA GeForce RTX 2080 Ti)

Iter  Time      Objective  L          Residual
-----------------------------------------------
   0  1.21e+00  7.795e+09  4.611e+03  4.126e+03
  10  2.33e+00  1.211e+06  4.611e+03  1.369e+01
  20  2.38e+00  2.102e+04  4.611e+03  1.476e+00
  30  2.43e+00  2.126e+03  4.611e+03  3.078e-01
  40  2.47e+00  8.819e+02  4.611e+03  9.760e-02
  50  2.52e+00  8.306e+02  4.611e+03  2.504e-02
  60  2.57e+00  8.240e+02  4.611e+03  1.229e-02
  70  2.62e+00  8.201e+02  4.611e+03  6.277e-03
  80  2.67e+00  8.197e+02  4.611e+03  2.978e-03
  90  2.72e+00  8.195e+02  4.611e+03  1.815e-03
  99  2.76e+00  8.195e+02  4.611e+03  1.059e-03

Plot the recovered coefficients and convergence statistics.

fig, ax = kplt.subplots(nrows=1, ncols=2, figsize=(12, 5))
kplt.plot(
    np.vstack((x_gt, x)).T,
    title="Coefficients",
    legend=("Ground Truth", "Recovered"),
    ax=ax[0],
)
kplt.plot(
    np.array((hist.Objective, hist.Residual)).T,
    ylog=True,
    title="Convergence",
    xlabel="Iteration",
    legend=("Objective", "Residual"),
    ax=ax[1],
)
fig.show()