.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/0005.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_0005.py: .. _gallery:0005: Cosine+Spikes, Dirac-Cosine Basis, Gaussian Measurements ============================================================= .. contents:: :depth: 2 :local: In this example we have #. A signal :math:`\by` consisting of a mixture of 3 cosine waves and 60 random spikes of total length 1024. #. A Dirac-Cosine two ortho basis :math:`\Psi` of shape 1024x2048. #. The sparse representation :math:`\bx` of the signal :math:`\by` in the basis :math:`\Psi` consisting of exactly 63 nonzero entries (corresponding to the spikes and the amplitudes of the cosine waves). #. A Gaussian sensing matrix :math:`\Phi` of shape 300x1024 making 300 random measurements in a vector :math:`\bb`. #. We are given :math:`\bb` and :math:`\bA = \Phi \Psi` and have to reconstruct :math:`\bx` using it. #. Then we can use :math:`\Psi` to compute :math:`\by = \Psi \bx`. .. math:: \bb = \bA \bx = \Phi \Psi \bx = \Phi \by. See also: * :ref:`api:problems` * :ref:`api:lop` .. GENERATED FROM PYTHON SOURCE LINES 35-45 .. 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.plots as crplot .. GENERATED FROM PYTHON SOURCE LINES 46-50 Setup ------------------------------ We shall construct our test signal and dictionary using our test problems module. .. GENERATED FROM PYTHON SOURCE LINES 50-56 .. code-block:: default from cr.sparse import problems prob = problems.generate('cosine-spikes:dirac-dct:gaussian', c=3, k=60) fig, ax = problems.plot(prob) .. image-sg:: /gallery/images/sphx_glr_0005_001.png :alt: Cosine with spikes, Cosine part of signal, Spikes part of signal, Dirac DCT Model, Measurements :srcset: /gallery/images/sphx_glr_0005_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 57-58 Let us access the relevant parts of our test problem .. GENERATED FROM PYTHON SOURCE LINES 58-68 .. code-block:: default # The combined linear operator (sensing matrix + dictionary) A = prob.A # The sparse representation of the signal in the dictionary x0 = prob.x # The Cosine+Spikes signal y0 = prob.y # The measurements b0 = prob.b .. GENERATED FROM PYTHON SOURCE LINES 69-71 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 71-73 .. code-block:: default print(crn.num_largest_coeffs_for_energy_percent(x0, 99.9)) .. rst-class:: sphx-glr-script-out .. code-block:: none 37 .. GENERATED FROM PYTHON SOURCE LINES 74-79 This number gives us an idea about the required sparsity to be configured for greedy pursuit algorithms. Although the exact sparsity of this representation is 63 but several of the spikes are too small and could be ignored for a reasonably good approximation. .. GENERATED FROM PYTHON SOURCE LINES 81-84 Sparse Recovery using Subspace Pursuit ------------------------------------------- We shall use subspace pursuit to reconstruct the signal. .. GENERATED FROM PYTHON SOURCE LINES 84-88 .. code-block:: default import cr.sparse.pursuit.sp as sp # We will first try to estimate a 50-sparse representation sol = sp.solve(A, b0, 50) .. GENERATED FROM PYTHON SOURCE LINES 89-90 This utility function helps us quickly analyze the quality of reconstruction .. GENERATED FROM PYTHON SOURCE LINES 90-92 .. code-block:: default problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 300, n: 2048 b_norm: original: 43.984 reconstruction: 43.982 SNR: 40.62 dB x_norm: original: 47.155 reconstruction: 47.147 SNR: 39.24 dB y_norm: original: 46.969 reconstruction: 46.986 SNR: 39.26 dB Sparsity: original: 37, reconstructed: 36, overlap: 36, ratio: 0.973 Iterations: 20 .. GENERATED FROM PYTHON SOURCE LINES 93-94 We will now try to estimate a 75-sparse representation .. GENERATED FROM PYTHON SOURCE LINES 94-96 .. code-block:: default sol = sp.solve(A, b0, 75) problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 300, n: 2048 b_norm: original: 43.984 reconstruction: 43.984 SNR: 82.42 dB x_norm: original: 47.155 reconstruction: 47.155 SNR: 79.39 dB y_norm: original: 46.969 reconstruction: 46.969 SNR: 79.58 dB Sparsity: original: 37, reconstructed: 37, overlap: 37, ratio: 1.000 Iterations: 10 .. GENERATED FROM PYTHON SOURCE LINES 97-99 Let us check if we correctly decoded all the nonzero entries in the sparse representation x .. GENERATED FROM PYTHON SOURCE LINES 99-101 .. code-block:: default problems.analyze_solution(prob, sol, perc=100) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 300, n: 2048 b_norm: original: 43.984 reconstruction: 43.984 SNR: 82.42 dB x_norm: original: 47.155 reconstruction: 47.155 SNR: 79.39 dB y_norm: original: 46.969 reconstruction: 46.969 SNR: 79.58 dB Sparsity: original: 63, reconstructed: 75, overlap: 62, ratio: 0.827 Iterations: 10 .. GENERATED FROM PYTHON SOURCE LINES 102-103 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 103-104 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 105-106 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 106-107 .. code-block:: default y = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 108-109 The estimated measurements .. GENERATED FROM PYTHON SOURCE LINES 109-111 .. code-block:: default b = A.times(x) .. GENERATED FROM PYTHON SOURCE LINES 112-113 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 113-122 .. code-block:: default def plot_representations(x0, x): ax = crplot.h_plots(2, height=2) ax[0].stem(x0, markerfmt='.') ax[0].set_title('Original representation') ax[1].stem(x, markerfmt='.') ax[1].set_title('Reconstructed representation') plot_representations(x0, x) .. image-sg:: /gallery/images/sphx_glr_0005_002.png :alt: Original representation, Reconstructed representation :srcset: /gallery/images/sphx_glr_0005_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 123-124 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 124-132 .. code-block:: default def plot_signals(y0, y): ax = crplot.h_plots(2, height=2) ax[0].plot(y0) ax[0].set_title('Original signal') ax[1].plot(y) ax[1].set_title('Reconstructed signal') plot_signals(y0, y) .. image-sg:: /gallery/images/sphx_glr_0005_003.png :alt: Original signal, Reconstructed signal :srcset: /gallery/images/sphx_glr_0005_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 133-134 Let us visualize the original and reconstructed measurements .. GENERATED FROM PYTHON SOURCE LINES 134-143 .. code-block:: default def plot_measurments(b0, b): ax = crplot.h_plots(2, height=2) ax[0].plot(b0) ax[0].set_title('Original measurements') ax[1].plot(b) ax[1].set_title('Reconstructed measurements') plot_measurments(b0, b) .. image-sg:: /gallery/images/sphx_glr_0005_004.png :alt: Original measurements, Reconstructed measurements :srcset: /gallery/images/sphx_glr_0005_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 144-147 Sparse Recovery using Compressive Sampling Matching Pursuit --------------------------------------------------------------- We shall now use compressive sampling matching pursuit to reconstruct the signal. .. GENERATED FROM PYTHON SOURCE LINES 147-152 .. code-block:: default import cr.sparse.pursuit.cosamp as cosamp # We will try to estimate a 75-sparse representation sol = cosamp.solve(A, b0, 75) problems.analyze_solution(prob, sol, perc=100) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 300, n: 2048 b_norm: original: 43.984 reconstruction: 43.984 SNR: 80.55 dB x_norm: original: 47.155 reconstruction: 47.155 SNR: 76.81 dB y_norm: original: 46.969 reconstruction: 46.969 SNR: 77.09 dB Sparsity: original: 63, reconstructed: 75, overlap: 62, ratio: 0.827 Iterations: 57 .. GENERATED FROM PYTHON SOURCE LINES 153-154 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 154-155 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 156-157 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 157-158 .. code-block:: default y = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 159-160 The estimated measurements .. GENERATED FROM PYTHON SOURCE LINES 160-163 .. code-block:: default b = A.times(x) .. GENERATED FROM PYTHON SOURCE LINES 164-165 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 165-168 .. code-block:: default plot_representations(x0, x) .. image-sg:: /gallery/images/sphx_glr_0005_005.png :alt: Original representation, Reconstructed representation :srcset: /gallery/images/sphx_glr_0005_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 169-170 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 170-172 .. code-block:: default plot_signals(y0, y) .. image-sg:: /gallery/images/sphx_glr_0005_006.png :alt: Original signal, Reconstructed signal :srcset: /gallery/images/sphx_glr_0005_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 173-174 Let us visualize the original and reconstructed measurements .. GENERATED FROM PYTHON SOURCE LINES 174-177 .. code-block:: default plot_measurments(b0, b) .. image-sg:: /gallery/images/sphx_glr_0005_007.png :alt: Original measurements, Reconstructed measurements :srcset: /gallery/images/sphx_glr_0005_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 178-180 Sparse Recovery using SPGL1 --------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 180-185 .. code-block:: default import cr.sparse.cvx.spgl1 as crspgl1 options = crspgl1.SPGL1Options(max_iters=1000) sol = crspgl1.solve_bp_jit(A, b0, options=options) problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 300, n: 2048 b_norm: original: 43.984 reconstruction: 43.984 SNR: 112.91 dB x_norm: original: 47.155 reconstruction: 46.907 SNR: 35.83 dB y_norm: original: 46.969 reconstruction: 46.781 SNR: 36.82 dB Sparsity: original: 37, reconstructed: 34, overlap: 34, ratio: 0.919 Iterations: 861 n_times: 1522, n_trans: 862 .. GENERATED FROM PYTHON SOURCE LINES 186-187 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 187-188 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 189-190 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 190-191 .. code-block:: default y = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 192-193 The estimated measurements .. GENERATED FROM PYTHON SOURCE LINES 193-196 .. code-block:: default b = A.times(x) .. GENERATED FROM PYTHON SOURCE LINES 197-198 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 198-201 .. code-block:: default plot_representations(x0, x) .. image-sg:: /gallery/images/sphx_glr_0005_008.png :alt: Original representation, Reconstructed representation :srcset: /gallery/images/sphx_glr_0005_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 202-203 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 203-205 .. code-block:: default plot_signals(y0, y) .. image-sg:: /gallery/images/sphx_glr_0005_009.png :alt: Original signal, Reconstructed signal :srcset: /gallery/images/sphx_glr_0005_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 206-207 Let us visualize the original and reconstructed measurements .. GENERATED FROM PYTHON SOURCE LINES 207-210 .. code-block:: default plot_measurments(b0, b) .. image-sg:: /gallery/images/sphx_glr_0005_010.png :alt: Original measurements, Reconstructed measurements :srcset: /gallery/images/sphx_glr_0005_010.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 211-225 Comments --------- * With K=50, SP recovery is slightly inaccurate. It also takes more (20) iterations to converge. * With K=75, SP is pretty good. It only missed one of the 63 nonzero entries. Also, SP converges in just 10 iterations. * With K=75, CoSaMP is also good but slightly poor. It also missed just one nonzero entry. But it seems like it missed a more significant entry compared to SP. Also, CoSaMP took 57 iterations to converge. * SPGL1 converges in 788 iterations to converge. Its model space and signal space SNR are not good. However, its measurement space SNR is pretty high. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 15.886 seconds) .. _sphx_glr_download_gallery_0005.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 0005.py <0005.py>` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 0005.ipynb <0005.ipynb>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_