Non-Negative Basis Pursuit DeNoising (APGM)

This example demonstrates the solution of a non-negative sparse coding problem

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

where \(D\) the dictionary, \(\mathbf{y}\) the signal to be represented, \(\mathbf{x}\) is the sparse representation, and \(\iota_{\mathrm{NN}}\) is the indicator function of the non-negativity constraint.

In this example the problem is solved via Accelerated PGM, using the proximal averaging method [62] to approximate the proximal operator of the sum of the \(\ell_1\) norm and an indicator function, while ADMM is used in a companion example.

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()

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

m = 32  # signal size
n = 128  # dictionary size
s = 10  # sparsity level

np.random.seed(1)
D = np.random.randn(m, n).astype(np.float32)
D = D / np.linalg.norm(D, axis=0, keepdims=True)  # normalize dictionary
L0 = max(np.linalg.norm(D, 2) ** 2, 5e1)

xt = np.zeros(n, dtype=np.float32)  # true signal
idx = np.random.randint(low=0, high=n, size=s)  # support of xt
xt[idx] = np.random.rand(s)
y = D @ xt + 5e-2 * np.random.randn(m)  # synthetic signal

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

Set up the forward operator and APGM solver object.

lmbda = 2e-1
A = linop.MatrixOperator(D)
f = loss.SquaredL2Loss(y=y, A=A)
g = functional.ProximalAverage([lmbda * functional.L1Norm(), functional.NonNegativeIndicator()])
maxiter = 250  # number of APGM iterations
solver = AcceleratedPGM(
    f=f, g=g, L0=L0, x0=A.adj(y), maxiter=maxiter, itstat_options={"display": True, "period": 20}
)

Run the solver.

print(f"Solving on {device_info()}\n")
x = solver.solve()

Solving on GPU (NVIDIA GeForce RTX 2080 Ti)

Iter  Time      Objective  L          Residual
-----------------------------------------------
   0  1.26e+00  1.835e+01  5.000e+01  1.558e+00
  20  4.75e+00  6.021e-01  5.000e+01  2.050e-02
  40  6.30e+00  4.900e-01  5.000e+01  2.644e-03
  60  7.84e+00  4.763e-01  5.000e+01  1.417e-03
  80  9.40e+00  4.755e-01  5.000e+01  1.393e-03
 100  1.08e+01  4.751e-01  5.000e+01  3.486e-04
 120  1.24e+01  4.749e-01  5.000e+01  1.473e-04
 140  1.37e+01  4.750e-01  5.000e+01  1.638e-04
 160  1.50e+01  4.748e-01  5.000e+01  8.783e-05
 180  1.65e+01  4.748e-01  5.000e+01  7.968e-05
 200  1.81e+01  4.748e-01  5.000e+01  5.151e-05
 220  1.96e+01  4.749e-01  5.000e+01  5.227e-05
 240  2.11e+01  4.749e-01  5.000e+01  3.555e-05
 249  2.20e+01  4.748e-01  5.000e+01  2.624e-05

Plot the recovered coefficients and signal.

fig, ax = kplt.subplots(nrows=1, ncols=2, figsize=(12, 5))
kplt.plot(
    np.vstack((xt, solver.x)).T,
    title="Coefficients",
    legend=("Ground Truth", "Recovered"),
    ax=ax[0],
)
kplt.plot(
    np.vstack((D @ xt, y, D @ solver.x)).T,
    title="Signal",
    legend=("Ground Truth", "Noisy", "Recovered"),
    ax=ax[1],
)
fig.show()