.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/0200_cs/sparse_binary_sensor.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_0200_cs_sparse_binary_sensor.py: .. _gallery:cs:sparse_binary_sensor: Sparse Binary Sensing Matrices =============================================== .. contents:: :depth: 2 :local: A (random) sparse binary sensing matrix has a very simple design. Assume that the signal space is :math:`\RR^N` and the measurement space is :math:`\RR^M`. Every column of a sparse binary sensing matrix has a 1 in exactly :math:`d` positions and 0s elsewhere. The indices at which ones are present are randomly selected for each column. Following is an example sparse binary matrix with 3 ones in each column .. math:: \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 0 & 1 & 1\\ 1 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0\\ 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0\\ 0 & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 1\\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ \end{bmatrix} From the perspective of algorithm design, we often require that the sensing matrix have unit norm columns. This can be easily attained for sparse binary matrices by scaling them with :math:`\frac{1}{\sqrt{d}}`. JAX provides an efficient way of storing sparse matrices in BCOO format. By default we employ this storage format for the (random) sparse binary matrices. .. GENERATED FROM PYTHON SOURCE LINES 50-51 Necessary imports .. GENERATED FROM PYTHON SOURCE LINES 51-63 .. code-block:: default import math import cr.nimble as crn import cr.sparse as crs import cr.sparse.dict as crdict import cr.sparse.data as crdata import cr.sparse.lop as crlop import cr.sparse.plots as crplots import numpy as np import jax import jax.numpy as jnp from jax import random .. GENERATED FROM PYTHON SOURCE LINES 64-65 Some random number generation keys .. GENERATED FROM PYTHON SOURCE LINES 65-68 .. code-block:: default key = random.PRNGKey(3) keys = random.split(key, 5) .. GENERATED FROM PYTHON SOURCE LINES 69-71 Creating Sparse Binary Sensing Matrices ---------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 71-82 .. code-block:: default # Matrix parameters m = 10 n = 16 d = 4 # construct a sparse binary sensing matrix A = crdict.sparse_binary_mtx(keys[0], m, n, d, normalize_atoms=False) # It is stored in a compressed BCOO format print(A) .. rst-class:: sphx-glr-script-out .. code-block:: none BCOO(uint8[10, 16], nse=64) .. GENERATED FROM PYTHON SOURCE LINES 83-84 If we wish to see its contents .. GENERATED FROM PYTHON SOURCE LINES 84-87 .. code-block:: default Ad = A.todense() print(Ad) .. rst-class:: sphx-glr-script-out .. code-block:: none [[0 0 0 1 1 0 0 1 1 0 0 1 1 1 1 0] [1 0 1 0 0 1 1 1 0 1 0 0 1 0 1 0] [0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1] [1 1 0 0 1 0 0 1 0 0 1 0 1 0 0 0] [1 1 1 0 1 1 0 0 1 1 1 0 0 0 0 0] [0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1] [0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0] [0 1 1 0 0 1 1 0 0 1 1 1 1 1 1 1] [0 0 0 1 1 1 1 1 1 0 0 0 0 1 0 1] [1 0 1 0 0 0 1 0 0 1 1 1 0 1 0 0]] .. GENERATED FROM PYTHON SOURCE LINES 88-89 We can quickly check that all columns have d ones .. GENERATED FROM PYTHON SOURCE LINES 89-92 .. code-block:: default print(jnp.sum(Ad, 0)) .. rst-class:: sphx-glr-script-out .. code-block:: none [4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4] .. GENERATED FROM PYTHON SOURCE LINES 93-96 By convention, we generate normalized sensing matrices by default. However, in the case of sparse binary matrices, it is more efficient to work with the unnormalized sensing matrix .. GENERATED FROM PYTHON SOURCE LINES 96-100 .. code-block:: default A = crdict.sparse_binary_mtx(keys[0], m, n, d) print(A.todense()) .. rst-class:: sphx-glr-script-out .. code-block:: none [[0. 0. 0. 0.5 0.5 0. 0. 0.5 0.5 0. 0. 0.5 0.5 0.5 0.5 0. ] [0.5 0. 0.5 0. 0. 0.5 0.5 0.5 0. 0.5 0. 0. 0.5 0. 0.5 0. ] [0. 0. 0. 0. 0. 0. 0. 0. 0.5 0. 0. 0.5 0. 0. 0. 0.5] [0.5 0.5 0. 0. 0.5 0. 0. 0.5 0. 0. 0.5 0. 0.5 0. 0. 0. ] [0.5 0.5 0.5 0. 0.5 0.5 0. 0. 0.5 0.5 0.5 0. 0. 0. 0. 0. ] [0. 0.5 0. 0.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.5] [0. 0. 0. 0.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.5 0. ] [0. 0.5 0.5 0. 0. 0.5 0.5 0. 0. 0.5 0.5 0.5 0.5 0.5 0.5 0.5] [0. 0. 0. 0.5 0.5 0.5 0.5 0.5 0.5 0. 0. 0. 0. 0.5 0. 0.5] [0.5 0. 0.5 0. 0. 0. 0.5 0. 0. 0.5 0.5 0.5 0. 0.5 0. 0. ]] .. GENERATED FROM PYTHON SOURCE LINES 101-106 Sparse Binary Sensing Linear Operators ---------------------------------------- It is often advantageous to work with matrices wrapped in our linear operator design. Let us construct the sensing matrix as a linear operator .. GENERATED FROM PYTHON SOURCE LINES 106-108 .. code-block:: default T = crlop.sparse_binary_dict(keys[0], m, n, d, normalize_atoms=False) .. GENERATED FROM PYTHON SOURCE LINES 109-110 We can extract the contents by multiplying with an identity matrix .. GENERATED FROM PYTHON SOURCE LINES 110-113 .. code-block:: default print(T.times(jnp.eye(n, dtype=int))) .. rst-class:: sphx-glr-script-out .. code-block:: none [[0 0 0 1 1 0 0 1 1 0 0 1 1 1 1 0] [1 0 1 0 0 1 1 1 0 1 0 0 1 0 1 0] [0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1] [1 1 0 0 1 0 0 1 0 0 1 0 1 0 0 0] [1 1 1 0 1 1 0 0 1 1 1 0 0 0 0 0] [0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1] [0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0] [0 1 1 0 0 1 1 0 0 1 1 1 1 1 1 1] [0 0 0 1 1 1 1 1 1 0 0 0 0 1 0 1] [1 0 1 0 0 0 1 0 0 1 1 1 0 1 0 0]] .. GENERATED FROM PYTHON SOURCE LINES 114-116 We can keep the normalization of the sensing matrix has a separate scaling operator .. GENERATED FROM PYTHON SOURCE LINES 116-125 .. code-block:: default # Let us construct a scaling operator d_scale = 1/ math.sqrt(d) T_scale = crlop.scalar_mult(d_scale, m) # Let us combine the scaling operator with the unnormalized sensing operator T_normed = crlop.compose(T_scale, T) # Verify the normalized operator print(T_normed.times(jnp.eye(n, dtype=int))) .. rst-class:: sphx-glr-script-out .. code-block:: none [[0. 0. 0. 0.5 0.5 0. 0. 0.5 0.5 0. 0. 0.5 0.5 0.5 0.5 0. ] [0.5 0. 0.5 0. 0. 0.5 0.5 0.5 0. 0.5 0. 0. 0.5 0. 0.5 0. ] [0. 0. 0. 0. 0. 0. 0. 0. 0.5 0. 0. 0.5 0. 0. 0. 0.5] [0.5 0.5 0. 0. 0.5 0. 0. 0.5 0. 0. 0.5 0. 0.5 0. 0. 0. ] [0.5 0.5 0.5 0. 0.5 0.5 0. 0. 0.5 0.5 0.5 0. 0. 0. 0. 0. ] [0. 0.5 0. 0.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.5] [0. 0. 0. 0.5 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.5 0. ] [0. 0.5 0.5 0. 0. 0.5 0.5 0. 0. 0.5 0.5 0.5 0.5 0.5 0.5 0.5] [0. 0. 0. 0.5 0.5 0.5 0.5 0.5 0.5 0. 0. 0. 0. 0.5 0. 0.5] [0.5 0. 0.5 0. 0. 0. 0.5 0. 0. 0.5 0.5 0.5 0. 0.5 0. 0. ]] .. GENERATED FROM PYTHON SOURCE LINES 126-130 Compressive Sensing --------------------- We shall use a larger problem to demonstrate the sensing capabilities of the sparse binary sensing operator. .. GENERATED FROM PYTHON SOURCE LINES 130-140 .. code-block:: default # Signal space dimension n = 1024 # Measurement space dimension m = 256 # Number of ones in each column d = 16 # Sparsity level of the signal to be sensed k = 40 .. GENERATED FROM PYTHON SOURCE LINES 141-144 Let us construct the unnormalized as well as normalized sensing operators. We shall use the unnormalized one during sensing but the normalized one during reconstruction. .. GENERATED FROM PYTHON SOURCE LINES 144-149 .. code-block:: default d_scale = 1/ math.sqrt(d) T = crlop.sparse_binary_dict(keys[1], m, n, d, normalize_atoms=False) T_scale = crlop.scalar_mult(d_scale, m) T_normed = crlop.compose(T_scale, T) .. GENERATED FROM PYTHON SOURCE LINES 150-151 We can quickly visualize the sparsity pattern of this sensing matrix .. GENERATED FROM PYTHON SOURCE LINES 151-155 .. code-block:: default A = T.times(jnp.eye(n, dtype=int)) ax = crplots.one_plot() ax.spy(A); .. image-sg:: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_001.png :alt: sparse binary sensor :srcset: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 156-157 Let us construct a sparse signal .. GENERATED FROM PYTHON SOURCE LINES 157-161 .. code-block:: default x, omega = crdata.sparse_normal_representations(keys[2], n, k) ax = crplots.one_plot() crplots.plot_sparse_signal(ax, x) .. image-sg:: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_002.png :alt: sparse binary sensor :srcset: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 162-167 .. code-block:: default y = T.times(x) ax = crplots.one_plot() crplots.plot_sparse_signal(ax, y) .. image-sg:: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_003.png :alt: sparse binary sensor :srcset: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 168-171 Sparse Recovery --------------------- We shall use various algorithms for reconstructing the original signal .. GENERATED FROM PYTHON SOURCE LINES 173-175 CoSaMP '''''''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 175-189 .. code-block:: default # Import the algorithm from cr.sparse.pursuit import cosamp # Solve the problem sol = cosamp.operator_solve_jit(T_normed, y, k) print(sol) # We need to scale the solution since the measurements were unscaled x_hat = sol.x * d_scale # Compute the SNR and PRD print(f'SNR: {crn.signal_noise_ratio(x, x_hat):.2f} dB, PRD: {crn.percent_rms_diff(x, x_hat):.0f} %') # Plot the original and the reconstructed signal ax = crplots.h_plots(2) crplots.plot_sparse_signals(ax, x, x_hat) .. image-sg:: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_004.png :alt: sparse binary sensor :srcset: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none iterations 7 m=256, n=1024, k=40 r_norm 4.468207e-14 x_norm 2.590201e+01 SNR: 296.00 dB, PRD: 0 % .. GENERATED FROM PYTHON SOURCE LINES 190-193 Subspace Pursuit '''''''''''''''''''''''''' Import the algorithm .. GENERATED FROM PYTHON SOURCE LINES 193-206 .. code-block:: default from cr.sparse.pursuit import sp # Solve the problem sol = sp.operator_solve_jit(T_normed, y, k) print(sol) # We need to scale the solution since the measurements were unscaled x_hat = sol.x * d_scale # Compute the SNR and PRD print(f'SNR: {crn.signal_noise_ratio(x, x_hat):.2f} dB, PRD: {crn.percent_rms_diff(x, x_hat):.0f} %') # Plot the original and the reconstructed signal ax = crplots.h_plots(2) crplots.plot_sparse_signals(ax, x, x_hat) .. image-sg:: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_005.png :alt: sparse binary sensor :srcset: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none iterations 6 m=256, n=1024, k=40 r_norm 4.307734e-14 x_norm 2.590201e+01 SNR: 295.15 dB, PRD: 0 % .. GENERATED FROM PYTHON SOURCE LINES 207-210 Hard Thresholding Pursuit '''''''''''''''''''''''''' Import the algorithm .. GENERATED FROM PYTHON SOURCE LINES 210-223 .. code-block:: default from cr.sparse.pursuit import htp # Solve the problem sol = htp.operator_solve_jit(T_normed, y, k) print(sol) # We need to scale the solution since the measurements were unscaled x_hat = sol.x * d_scale # Compute the SNR and PRD print(f'SNR: {crn.signal_noise_ratio(x, x_hat):.2f} dB, PRD: {crn.percent_rms_diff(x, x_hat):.0f} %') # Plot the original and the reconstructed signal ax = crplots.h_plots(2) crplots.plot_sparse_signals(ax, x, x_hat) .. image-sg:: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_006.png :alt: sparse binary sensor :srcset: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none iterations 20 m=256, n=1024, k=40 r_norm 2.904729e-02 x_norm 2.590449e+01 SNR: 56.77 dB, PRD: 0 % .. GENERATED FROM PYTHON SOURCE LINES 224-227 Truncated Newton Interior Points Method '''''''''''''''''''''''''''''''''''''''''''''''' Import the algorithm .. GENERATED FROM PYTHON SOURCE LINES 227-239 .. code-block:: default from cr.sparse.cvx import l1ls # Solve the problem # Note that this algorithm doesn't require sparsity level k as input sol = l1ls.solve_jit(T_normed, y, 1e-2) print(sol) # We need to scale the solution since the measurements were unscaled x_hat = sol.x * d_scale # Compute the SNR and PRD print(f'SNR: {crn.signal_noise_ratio(x, x_hat):.2f} dB, PRD: {crn.percent_rms_diff(x, x_hat):.0f} %') # Plot the original and the reconstructed signal ax = crplots.h_plots(2) crplots.plot_sparse_signals(ax, x, x_hat) .. image-sg:: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_007.png :alt: sparse binary sensor :srcset: /gallery/0200_cs/images/sphx_glr_sparse_binary_sensor_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none iterations 19 n_times 1121 n_trans 1121 r_norm 3.622522e-02 SNR: 51.38 dB, PRD: 0 % .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 13.591 seconds) .. _sphx_glr_download_gallery_0200_cs_sparse_binary_sensor.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: sparse_binary_sensor.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: sparse_binary_sensor.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_