.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/0002.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_0002.py: .. _gallery:0002: Blocks Signal in Haar Basis ============================= .. contents:: :depth: 2 :local: A Blocks signal as proposed by Donoho et al. in Wavelab :cite:`buckheit1995wavelab` is a concatenation of blocks with different heights. It has a sparse representation in a Haar basis. In this test problem with the structure :math:`\bb = \bA \bx`, the signal :math:`\bb` is the blocks signal (of length 1024), :math:`\bA` is a Haar basis with 5 levels of decomposition and :math:`\bx` is the sparse representation of the signal in this basis. This basis is real, complete and orthonormal. Hence, a simple solution for getting :math:`\bx` from :math:`\bb` is the formula: .. math:: \bx = \bA^T \bb. This test problem is useful in identifying basic mistakes in a sparse recovery algorithm. This problem should be very easy to solve by any sparse recovery algorithm. However, there is a caveat. If a sparse recovery algorithm depends on the expected sparsity of the signal (typically a parameter K in greedy algorithms), the reconstruction would fail if K is specified below the actual number of significant components of :math:`\bx`. See also: * :ref:`api:problems` * :ref:`api:lop` .. GENERATED FROM PYTHON SOURCE LINES 47-58 .. code-block:: default # Configure JAX to work with 64-bit floating point precision. from jax.config import config config.update("jax_enable_x64", True) import jax.numpy as jnp import cr.nimble as crn import cr.sparse as crs import cr.sparse.plots as crplot .. GENERATED FROM PYTHON SOURCE LINES 59-63 Setup ------------------------------ We shall construct our test signal and dictionary using our test problems module. .. GENERATED FROM PYTHON SOURCE LINES 63-69 .. code-block:: default from cr.sparse import problems prob = problems.generate('blocks:haar') fig, ax = problems.plot(prob) .. image-sg:: /gallery/images/sphx_glr_0002_001.png :alt: Block signal, Haar coefficients of the blocks signal :srcset: /gallery/images/sphx_glr_0002_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 70-71 Let us access the relevant parts of our test problem .. GENERATED FROM PYTHON SOURCE LINES 71-80 .. code-block:: default # The sparsifying basis linear operator A = prob.A # The Blocks signal b0 = prob.b # The sparse representation of the Blocks signal in the dictionary x0 = prob.x .. GENERATED FROM PYTHON SOURCE LINES 81-83 Check how many coefficients in the sparse representation are sufficient to capture 99.9% of the energy of the signal .. GENERATED FROM PYTHON SOURCE LINES 83-85 .. code-block:: default print(crn.num_largest_coeffs_for_energy_percent(x0, 99.9)) .. rst-class:: sphx-glr-script-out .. code-block:: none 64 .. GENERATED FROM PYTHON SOURCE LINES 86-88 This number gives us an idea about the required sparsity to be configured for greedy pursuit algorithms. .. GENERATED FROM PYTHON SOURCE LINES 90-93 Sparse Recovery using Subspace Pursuit ------------------------------------------- We shall use subspace pursuit to reconstruct the signal. .. GENERATED FROM PYTHON SOURCE LINES 93-96 .. code-block:: default import cr.sparse.pursuit.sp as sp # We will try to estimate a 100-sparse representation sol = sp.solve(A, b0, 100) .. GENERATED FROM PYTHON SOURCE LINES 97-98 This utility function helps us quickly analyze the quality of reconstruction .. GENERATED FROM PYTHON SOURCE LINES 98-100 .. code-block:: default problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 1024 b_norm: original: 78.899 reconstruction: 78.899 SNR: 289.74 dB x_norm: original: 78.899 reconstruction: 78.899 SNR: 135.52 dB Sparsity: original: 64, reconstructed: 64, overlap: 64, ratio: 1.000 Iterations: 1 .. GENERATED FROM PYTHON SOURCE LINES 101-102 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 102-103 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 104-105 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 105-107 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 108-109 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 109-115 .. code-block:: default ax = crplot.h_plots(2) ax[0].stem(x0, markerfmt='.') ax[1].stem(x, markerfmt='.') .. image-sg:: /gallery/images/sphx_glr_0002_002.png :alt: 0002 :srcset: /gallery/images/sphx_glr_0002_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 116-117 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 117-123 .. code-block:: default ax = crplot.h_plots(2) ax[0].plot(b0) ax[1].plot(b) .. image-sg:: /gallery/images/sphx_glr_0002_003.png :alt: 0002 :srcset: /gallery/images/sphx_glr_0002_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 124-127 Sparse Recovery using Compressive Sampling Matching Pursuit --------------------------------------------------------------- We shall now use compressive sampling matching pursuit to reconstruct the signal. .. GENERATED FROM PYTHON SOURCE LINES 127-132 .. code-block:: default import cr.sparse.pursuit.cosamp as cosamp # We will try to estimate a 100-sparse representation sol = cosamp.solve(A, b0, 100) problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 1024 b_norm: original: 78.899 reconstruction: 78.899 SNR: 289.03 dB x_norm: original: 78.899 reconstruction: 78.899 SNR: 135.52 dB Sparsity: original: 64, reconstructed: 64, overlap: 64, ratio: 1.000 Iterations: 2 .. GENERATED FROM PYTHON SOURCE LINES 133-134 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 134-135 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 136-137 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 137-140 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 141-142 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 142-147 .. code-block:: default ax = crplot.h_plots(2) ax[0].stem(x0, markerfmt='.') ax[1].stem(x, markerfmt='.') .. image-sg:: /gallery/images/sphx_glr_0002_004.png :alt: 0002 :srcset: /gallery/images/sphx_glr_0002_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 148-149 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 149-154 .. code-block:: default ax = crplot.h_plots(2) ax[0].plot(b0) ax[1].plot(b) .. image-sg:: /gallery/images/sphx_glr_0002_005.png :alt: 0002 :srcset: /gallery/images/sphx_glr_0002_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 155-157 Sparse Recovery using SPGL1 --------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 157-161 .. code-block:: default import cr.sparse.cvx.spgl1 as crspgl1 options = crspgl1.SPGL1Options() tracker = crs.ProgressTracker(x0=x0) sol = crspgl1.solve_bp_jit(A, b0, options=options, tracker=tracker) .. rst-class:: sphx-glr-script-out .. code-block:: none [1] x_norm: 0.00e+00, r_norm: 7.89e+01, SNR: 0.00 dB [2] x_norm: 3.05e+01, r_norm: 4.84e+01, SNR: 4.25 dB [3] x_norm: 5.74e+01, r_norm: 3.08e+01, SNR: 8.17 dB [4] x_norm: 5.74e+01, r_norm: 3.08e+01, SNR: 8.17 dB [5] x_norm: 7.61e+01, r_norm: 4.28e+00, SNR: 25.30 dB [6] x_norm: 7.61e+01, r_norm: 4.28e+00, SNR: 25.30 dB [7] x_norm: 7.89e+01, r_norm: 2.48e-07, SNR: 135.63 dB Algorithm converged in 7 iterations. .. GENERATED FROM PYTHON SOURCE LINES 162-163 Let's analyze the solution quality .. GENERATED FROM PYTHON SOURCE LINES 163-165 .. code-block:: default problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 1024 b_norm: original: 78.899 reconstruction: 78.899 SNR: 170.06 dB x_norm: original: 78.899 reconstruction: 78.899 SNR: 135.63 dB Sparsity: original: 64, reconstructed: 64, overlap: 64, ratio: 1.000 Iterations: 7 n_times: 10, n_trans: 8 .. GENERATED FROM PYTHON SOURCE LINES 166-167 Let's plot the progress of SPGL1 over different iterations .. GENERATED FROM PYTHON SOURCE LINES 167-170 .. code-block:: default ax = crplot.one_plot(height=6) tracker.plot_progress(ax) .. image-sg:: /gallery/images/sphx_glr_0002_006.png :alt: 0002 :srcset: /gallery/images/sphx_glr_0002_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 171-172 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 172-173 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 174-175 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 175-178 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 179-180 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 180-185 .. code-block:: default ax = crplot.h_plots(2) ax[0].stem(x0, markerfmt='.') ax[1].stem(x, markerfmt='.') .. image-sg:: /gallery/images/sphx_glr_0002_007.png :alt: 0002 :srcset: /gallery/images/sphx_glr_0002_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 186-187 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 187-192 .. code-block:: default ax = crplot.h_plots(2) ax[0].plot(b0) ax[1].plot(b) .. image-sg:: /gallery/images/sphx_glr_0002_008.png :alt: 0002 :srcset: /gallery/images/sphx_glr_0002_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 193-197 Comments --------- We note that SP converges in a single iteration, CoSaMP takes two iterations, while SPGL1 takes 9 iterations to converge. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 15.392 seconds) .. _sphx_glr_download_gallery_0002.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 0002.py <0002.py>` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 0002.ipynb <0002.ipynb>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_