.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/0300_cvx/spikes_l1ls.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_0300_cvx_spikes_l1ls.py: Recovering spikes via TNIPM ===================================== This example has following features: * A sparse signal consists of a small number of spikes. * The sensing matrix is a random dictionary with orthonormal rows. * The number of measurements is one fourth of ambient dimensions. * The measurements are corrupted by noise. * Truncated Newton Interior Points Method (TNIPM) a.k.a. l1-ls algorithm is being used for recovery. .. GENERATED FROM PYTHON SOURCE LINES 18-19 Let's import necessary libraries .. GENERATED FROM PYTHON SOURCE LINES 19-36 .. code-block:: default import jax.numpy as jnp from jax import random norm = jnp.linalg.norm import matplotlib as mpl import matplotlib.pyplot as plt import cr.sparse as crs import cr.sparse.data as crdata import cr.sparse.lop as lop import cr.sparse.cvx.l1ls as l1ls from cr.nimble.dsp import ( hard_threshold_by, support, largest_indices_by ) .. GENERATED FROM PYTHON SOURCE LINES 37-39 Setup ------ .. GENERATED FROM PYTHON SOURCE LINES 39-51 .. code-block:: default # Number of measurements m = 2**10 # Ambient dimension n = 2**12 # Number of spikes (sparsity) k = 160 print(f'{m=}, {n=}') key = random.PRNGKey(0) keys = random.split(key, 4) .. rst-class:: sphx-glr-script-out .. code-block:: none m=1024, n=4096 .. GENERATED FROM PYTHON SOURCE LINES 52-54 The Spikes -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 54-58 .. code-block:: default xs, omega = crdata.sparse_spikes(keys[0], n, k) plt.figure(figsize=(8,6), dpi= 100, facecolor='w', edgecolor='k') plt.plot(xs) .. image-sg:: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_001.png :alt: spikes l1ls :srcset: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 59-61 The Sparsifying Basis -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 61-63 .. code-block:: default A = lop.random_orthonormal_rows_dict(keys[1], m, n) .. GENERATED FROM PYTHON SOURCE LINES 64-66 Measurement process ------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 66-77 .. code-block:: default # Clean measurements bs = A.times(xs) # Noise sigma = 0.01 noise = sigma * random.normal(keys[2], (m,)) # Noisy measurements b = bs + noise plt.figure(figsize=(8,6), dpi= 100, facecolor='w', edgecolor='k') plt.plot(b) .. image-sg:: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_002.png :alt: spikes l1ls :srcset: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 78-80 Recovery using TNIPM ------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 80-104 .. code-block:: default # We need to estimate the regularization parameter Atb = A.trans(b) tau = float(0.1 * jnp.max(jnp.abs(Atb))) print(f'{tau=}') # Now run the solver sol = l1ls.solve_jit(A, b, tau) # number of L1-LS iterations iterations = int(sol.iterations) # number of A x operations n_times = int(sol.n_times) # number of A^H y operations n_trans = int(sol.n_trans) print(f'{iterations=} {n_times=} {n_trans=}') # residual norm r_norm = norm(sol.x) print(f'{r_norm=:.3e}') # relative error rel_error = norm(xs - sol.x) / norm(xs) print(f'{rel_error=:.3e}') .. rst-class:: sphx-glr-script-out .. code-block:: none tau=0.05101185632358787 iterations=17 n_times=173 n_trans=174 r_norm=1.122e+01 rel_error=1.224e-01 .. GENERATED FROM PYTHON SOURCE LINES 105-107 Solution ------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 107-113 .. code-block:: default plt.figure(figsize=(8,6), dpi= 100, facecolor='w', edgecolor='k') plt.subplot(211) plt.plot(xs) plt.subplot(212) plt.plot(sol.x) .. image-sg:: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_003.png :alt: spikes l1ls :srcset: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 114-116 The magnitudes of non-zero values ''''''''''''''''''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 116-119 .. code-block:: default plt.figure(figsize=(8,6), dpi= 100, facecolor='w', edgecolor='k') plt.plot(jnp.sort(jnp.abs(sol.x))) .. image-sg:: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_004.png :alt: spikes l1ls :srcset: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 120-122 Thresholding for large values ''''''''''''''''''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 122-129 .. code-block:: default x = hard_threshold_by(sol.x, 0.5) plt.figure(figsize=(8,6), dpi= 100, facecolor='w', edgecolor='k') plt.subplot(211) plt.plot(xs) plt.subplot(212) plt.plot(x) .. image-sg:: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_005.png :alt: spikes l1ls :srcset: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 130-132 Verifying the support recovery ''''''''''''''''''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 132-137 .. code-block:: default support_xs = support(xs) support_x = support(x) jnp.all(jnp.equal(support_xs, support_x)) .. rst-class:: sphx-glr-script-out .. code-block:: none Array(True, dtype=bool) .. GENERATED FROM PYTHON SOURCE LINES 138-140 Improvement using least squares over support ------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 140-161 .. code-block:: default # Identify the sub-matrix of columns for the support of recovered solution's large entries support_x = largest_indices_by(sol.x, 0.5) AI = A.columns(support_x) print(AI.shape) # Solve the least squares problem over these columns x_I, residuals, rank, s = jnp.linalg.lstsq(AI, b) # fill the non-zero entries into the sparse least squares solution x_ls = jnp.zeros_like(xs) x_ls = x_ls.at[support_x].set(x_I) # relative error ls_rel_error = norm(xs - x_ls) / norm(xs) print(f'{ls_rel_error=:.3e}') plt.figure(figsize=(8,6), dpi= 100, facecolor='w', edgecolor='k') plt.subplot(211) plt.plot(xs) plt.subplot(212) plt.plot(x_ls) .. image-sg:: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_006.png :alt: spikes l1ls :srcset: /gallery/0300_cvx/images/sphx_glr_spikes_l1ls_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (1024, 160) ls_rel_error=2.070e-02 [] .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 7.767 seconds) .. _sphx_glr_download_gallery_0300_cvx_spikes_l1ls.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: spikes_l1ls.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: spikes_l1ls.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_