.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/0001.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_0001.py: .. _gallery:0001: HeaviSine Signal in Fourier-HeaviSide Basis ============================================ .. contents:: :depth: 2 :local: A HeaviSine signal as proposed by Donoho et al. in Wavelab :cite:`buckheit1995wavelab` is a sign wave with jump discontinuities. It is not sparse in standard (dirac) basis. It also doesn't have a sparse representation in the Fourier basis. However, it does have a sparse representation in the Fourier-Heaviside dictionary. In this example, we will generate a test problem of HeaviSine signal and construct its sparse representation through various algorithms. The dictionary is a concatenation of the (orthonormal) Fourier basis and a (non-orthogonal) HeaviSide basis. HeaviSide basis is derived from the `HeaviSide step function `_. In its finite dimensional discrete form it looks like an :math:`n \times n` matrix that has ones below and on the diagonal and zeros elsewhere. In other words, all the elements above the diagonal are zero and rest are one. Typical sparse reconstruction algorithms assume that the atoms in a sparsifying dictionary have unit norms. We provide both the unnormalized (with one in lower triangular part) and normalized versions of the HeaviSide basis in ``cr.sparse.lop`` module. In this example, we shall use the normalized HeaviSide basis. Since the dictionary contains a Fourier basis, hence the representation of the HeaviSine signal in this dictionary is a complex valued representation. The signal itself however is real. See also: * :ref:`api:problems` * :ref:`api:lop` * :ref:`api:l1min` .. GENERATED FROM PYTHON SOURCE LINES 54-62 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 63-67 Setup ------------------------------ We shall construct our test signal and dictionary using our test problems module. .. GENERATED FROM PYTHON SOURCE LINES 67-72 .. code-block:: default from cr.sparse import problems prob = problems.generate('heavi-sine:fourier:heavi-side') fig, ax = problems.plot(prob) .. image-sg:: /gallery/images/sphx_glr_0001_001.png :alt: HeaviSine signal, Sine part of signal, Piecewise constant component of the HeaviSine signal, Real part of Fourier coefficients for the Sine component, Imaginary part of Fourier coefficients for the Sine component, Coefficients for the jumps part in Heaviside (non-orthogonal) basis :srcset: /gallery/images/sphx_glr_0001_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 73-74 Let us access the relevant parts of our test problem .. GENERATED FROM PYTHON SOURCE LINES 74-83 .. code-block:: default # The sparsifying basis linear operator A = prob.A # The HeaviSine signal b0 = prob.b # The sparse representation of the HeaviSide signal in the dictionary x0 = prob.x .. GENERATED FROM PYTHON SOURCE LINES 84-86 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 86-88 .. code-block:: default print(crn.num_largest_coeffs_for_energy_percent(x0, 99.9)) .. rst-class:: sphx-glr-script-out .. code-block:: none 5 .. GENERATED FROM PYTHON SOURCE LINES 89-92 Sparse Recovery using Subspace Pursuit ------------------------------------------- We shall use subspace pursuit to reconstruct the signal. .. GENERATED FROM PYTHON SOURCE LINES 92-97 .. code-block:: default import cr.sparse.pursuit.sp as sp # We will try to estimate a 10-sparse representation sol = sp.solve(prob.A, prob.b, 10) print(sol) .. rst-class:: sphx-glr-script-out .. code-block:: none iterations 4 m=1024, n=2048, k=10 r_norm 9.102335e-14 x_norm 1.193985e+02 .. GENERATED FROM PYTHON SOURCE LINES 98-99 The sparse representation estimated by subspace pursuit .. GENERATED FROM PYTHON SOURCE LINES 99-100 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 101-102 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 102-104 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 105-106 Check if we could reconstruct the sparse representation correctly .. GENERATED FROM PYTHON SOURCE LINES 106-107 .. code-block:: default print(f'Model Space SNR: {crn.signal_noise_ratio(x0, x)} dB') .. rst-class:: sphx-glr-script-out .. code-block:: none Model Space SNR: 9.009642143585591 dB .. GENERATED FROM PYTHON SOURCE LINES 108-115 The SNR between the expected sparse representation and the recovered sparse representation is low. HeaviSide basis is highly correlated (high coherence). Hence, we couldn't make the exact recovery of the original sparse representation. Nevertheless we indeed recovered a good sparse representation since the residual norm is very small .. GENERATED FROM PYTHON SOURCE LINES 117-118 Check if we could reconstruct the signal correctly .. GENERATED FROM PYTHON SOURCE LINES 118-119 .. code-block:: default print(f'Signal Space SNR: {crn.signal_noise_ratio(b0, b)} dB') .. rst-class:: sphx-glr-script-out .. code-block:: none Signal Space SNR: 302.32500793997445 dB .. GENERATED FROM PYTHON SOURCE LINES 120-121 The reconstruction is indeed excellent. .. GENERATED FROM PYTHON SOURCE LINES 125-126 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 126-131 .. code-block:: default import cr.sparse.plots as crplot ax = crplot.h_plots(2) ax[0].plot(b0) ax[1].plot(b) .. image-sg:: /gallery/images/sphx_glr_0001_002.png :alt: 0001 :srcset: /gallery/images/sphx_glr_0001_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 132-134 Sparse Recovery using SPGL1 --------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 134-139 .. code-block:: default import cr.sparse.cvx.spgl1 as crspgl1 options = crspgl1.SPGL1Options(max_iters=2000) sol = crspgl1.solve_bp_jit(A, b0, options=options) problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 75.465 reconstruction: 75.465 SNR: 117.56 dB x_norm: original: 127.687 reconstruction: 71.466 SNR: 2.33 dB Sparsity: original: 5, reconstructed: 282, overlap: 4, ratio: 0.014 Iterations: 1508 n_times: 2636, n_trans: 1509 .. GENERATED FROM PYTHON SOURCE LINES 140-141 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 141-142 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 143-144 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 144-146 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 147-148 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 148-152 .. code-block:: default ax = crplot.h_plots(2) ax[0].plot(b0) ax[1].plot(b) .. image-sg:: /gallery/images/sphx_glr_0001_003.png :alt: 0001 :srcset: /gallery/images/sphx_glr_0001_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 153-159 Comments --------- * Subspace Pursuit converges very fast in 4 iterations. * SPGL1 takes 1600+ iterations to converge. Still it is inaccurate. * The HeaviSide basis is highly coherent. It doesn't satisfy RIP guarantees. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 10.706 seconds) .. _sphx_glr_download_gallery_0001.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 0001.py <0001.py>` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 0001.ipynb <0001.ipynb>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_