Compressive Sensing Operators

Let’s import necessary libraries

import jax.numpy as jnp
# Random numbers
from jax import random
# For plotting diagrams
import matplotlib.pyplot as plt
## CR-Sparse modules
import cr.nimble as cnb
# Linear operators
from cr.sparse import lop
# Some random number seeds
key = random.PRNGKey(0)
keys = random.split(key, 10)

Rademacher sensing operators

Size of the sensing matrix

m, n = 4, 8

Unnormalized atoms / columns

Construct the operator wrapping the sensing matrix

Phi = lop.rademacher_dict(keys[0], m, n, normalize_atoms=False)

Let’s print the contents of the matrix

print(lop.to_matrix(Phi))
[[-1.  1. -1. -1.  1. -1. -1.  1.]
 [ 1. -1. -1. -1.  1. -1.  1.  1.]
 [-1. -1. -1. -1. -1.  1. -1. -1.]
 [ 1. -1. -1.  1.  1. -1.  1. -1.]]

Let’s print the contents of the adjoint matrix

print(lop.to_adjoint_matrix(Phi))
[[-1.  1. -1.  1.]
 [ 1. -1. -1. -1.]
 [-1. -1. -1. -1.]
 [-1. -1. -1.  1.]
 [ 1.  1. -1.  1.]
 [-1. -1.  1. -1.]
 [-1.  1. -1.  1.]
 [ 1.  1. -1. -1.]]

Let’s apply the forward operator

x = jnp.array([3, -2, 4, 8, 2, 2, 6, -1])
y = Phi.times(x)
print(y)
[-24  -2 -18  16]

Let’s apply the adjoint operator

z = Phi.trans(y)
print(z)
[ 56 -20  28  60   8  -8  56 -24]

Let’s JIT compile the trans and times functions

Phi = lop.jit(Phi)

Let’s verify the forward operator after JIT compile

x = jnp.array([3, -2, 4, 8, 2, 2, 6, -1])
y = Phi.times(x)
print(y)
[-24  -2 -18  16]

Let’s verify the adjoint operator after JIT compile

z = Phi.trans(y)
print(z)
[ 56 -20  28  60   8  -8  56 -24]

Normalized atoms / columns

Construct the operator wrapping the sensing matrix

Phi = lop.rademacher_dict(keys[0], m, n, normalize_atoms=True)

Let’s print the contents of the matrix

print(lop.to_matrix(Phi))
[[-0.5  0.5 -0.5 -0.5  0.5 -0.5 -0.5  0.5]
 [ 0.5 -0.5 -0.5 -0.5  0.5 -0.5  0.5  0.5]
 [-0.5 -0.5 -0.5 -0.5 -0.5  0.5 -0.5 -0.5]
 [ 0.5 -0.5 -0.5  0.5  0.5 -0.5  0.5 -0.5]]

Check the column wise norms

print(cnb.norms_l2_cw(lop.to_matrix(Phi)))
[1. 1. 1. 1. 1. 1. 1. 1.]

Let’s print the contents of the adjoint matrix

print(lop.to_adjoint_matrix(Phi))
[[-0.5  0.5 -0.5  0.5]
 [ 0.5 -0.5 -0.5 -0.5]
 [-0.5 -0.5 -0.5 -0.5]
 [-0.5 -0.5 -0.5  0.5]
 [ 0.5  0.5 -0.5  0.5]
 [-0.5 -0.5  0.5 -0.5]
 [-0.5  0.5 -0.5  0.5]
 [ 0.5  0.5 -0.5 -0.5]]

Let’s apply the forward operator

x = jnp.array([3, -2, 4, 8, 2, 2, 6, -1])
y = Phi.times(x)
print(y)
[-12.  -1.  -9.   8.]

Let’s apply the adjoint operator

z = Phi.trans(y)
print(z)
[14. -5.  7. 15.  2. -2. 14. -6.]

Let’s JIT compile the trans and times functions

Phi = lop.jit(Phi)

Let’s verify the forward operator after JIT compile

y = Phi.times(x)
print(y)
[-12.  -1.  -9.   8.]

Let’s verify the adjoint operator after JIT compile

z = Phi.trans(y)
print(z)
[14. -5.  7. 15.  2. -2. 14. -6.]

Column wise application on input

Number of signals

k = 5

Prepare a signal matrix

X = random.randint(keys[1], (n, k), -5, 5)
print(X)
[[-4 -5 -2 -4  2]
 [ 0  2 -2  1  3]
 [ 2 -2  4 -4 -2]
 [ 2  4 -5  1 -1]
 [ 2 -5 -2 -2 -4]
 [-4  2 -3 -1 -4]
 [-1 -1 -1  4  1]
 [-3 -4 -1 -4  1]]

Let’s apply the forward operator column wise

Y = Phi.times(X)
print(Y)
[[  2.   -2.5   1.   -0.5   2. ]
 [ -3.  -10.5   0.   -1.5   2. ]
 [ -1.    6.5   3.    3.5  -2. ]
 [  2.   -2.5  -4.    3.5   0. ]]

Let’s apply the adjoint operator column wise

Z = Phi.trans(Y)
print(Z)
[[ -1.   -8.5  -4.   -0.5   1. ]
 [  2.    2.    1.   -3.    1. ]
 [  0.    4.5   0.   -2.5  -1. ]
 [  2.    2.   -4.    1.   -1. ]
 [  1.  -11.   -3.   -1.    3. ]
 [ -1.   11.    3.    1.   -3. ]
 [ -1.   -8.5  -4.   -0.5   1. ]
 [ -1.   -8.5   1.   -4.5   3. ]]

Row wise application on input

Construct the operator wrapping the sensing matrix

Phi = lop.rademacher_dict(keys[0], m, n, normalize_atoms=True, axis=1)

Prepare a signal matrix with each signal along a row

XT = X.T

Let’s apply the forward operator column wise

Y = Phi.times(XT)
print(Y)
[[  2.   -3.   -1.    2. ]
 [ -2.5 -10.5   6.5  -2.5]
 [  1.    0.    3.   -4. ]
 [ -0.5  -1.5   3.5   3.5]
 [  2.    2.   -2.    0. ]]

Let’s apply the adjoint operator column wise

Z = Phi.trans(Y)
print(Z)
[[ -1.    2.    0.    2.    1.   -1.   -1.   -1. ]
 [ -8.5   2.    4.5   2.  -11.   11.   -8.5  -8.5]
 [ -4.    1.    0.   -4.   -3.    3.   -4.    1. ]
 [ -0.5  -3.   -2.5   1.   -1.    1.   -0.5  -4.5]
 [  1.    1.   -1.   -1.    3.   -3.    1.    3. ]]

Square sensing matrix

Construct the operator wrapping the sensing matrix

Phi = lop.rademacher_dict(keys[0], n, normalize_atoms=False)

Let’s print the contents of the matrix

print(lop.to_matrix(Phi))
[[ 1.  1. -1. -1. -1. -1.  1.  1.]
 [-1.  1.  1.  1. -1. -1. -1. -1.]
 [-1. -1.  1.  1.  1. -1.  1. -1.]
 [ 1. -1.  1. -1.  1.  1.  1. -1.]
 [-1.  1. -1.  1.  1.  1. -1. -1.]
 [-1.  1. -1.  1.  1. -1.  1.  1.]
 [ 1. -1.  1. -1. -1. -1. -1.  1.]
 [-1. -1.  1. -1.  1. -1.  1. -1.]]

Gaussian sensing operators

Unnormalized atoms / columns

Construct the operator wrapping the sensing matrix

Phi = lop.gaussian_dict(keys[0], m, n, normalize_atoms=False)

Let’s print the contents of the matrix

print(lop.to_matrix(Phi))
[[-0.76261741 -0.54906996  0.35375743  0.44785326  0.20306852  0.4002937
  -0.10082346 -1.38919403]
 [ 0.68214811 -0.73785154 -0.31989581 -0.37379108  0.02430423  0.60747778
   0.24837508  0.39562076]
 [ 0.08833386  0.49151262 -0.45341975 -0.63881669 -0.23097626  0.10892717
  -0.75912382  0.09899961]
 [-0.97121466  0.38958316 -0.38142231  1.05170598 -0.59028129 -0.18775354
  -0.78933777  0.06830431]]

Let’s print the contents of the adjoint matrix

print(lop.to_adjoint_matrix(Phi))
[[-0.76261741  0.68214811  0.08833386 -0.97121466]
 [-0.54906996 -0.73785154  0.49151262  0.38958316]
 [ 0.35375743 -0.31989581 -0.45341975 -0.38142231]
 [ 0.44785326 -0.37379108 -0.63881669  1.05170598]
 [ 0.20306852  0.02430423 -0.23097626 -0.59028129]
 [ 0.4002937   0.60747778  0.10892717 -0.18775354]
 [-0.10082346  0.24837508 -0.75912382 -0.78933777]
 [-1.38919403  0.39562076  0.09899961  0.06830431]]

Let’s apply the forward operator

x = jnp.array([3, -2, 4, 8, 2, 2, 6, -1])
y = Phi.times(x)
print(y)
[  5.79912122   1.61042927 -12.54007687  -3.16525231]

Let’s apply the adjoint operator

z = Phi.trans(y)
print(z)
[ -1.3575335  -11.76911601   8.42952898   6.67708692   5.98160846
   2.52798389  11.83322722  -8.87664847]

Let’s JIT compile the trans and times functions

Phi = lop.jit(Phi)

Let’s verify the forward operator after JIT compile

x = jnp.array([3, -2, 4, 8, 2, 2, 6, -1])
y = Phi.times(x)
print(y)
[  5.79912122   1.61042927 -12.54007687  -3.16525231]

Let’s verify the adjoint operator after JIT compile

z = Phi.trans(y)
print(z)
[ -1.3575335  -11.76911601   8.42952898   6.67708692   5.98160846
   2.52798389  11.83322722  -8.87664847]

Normalized atoms / columns

Construct the operator wrapping the sensing matrix

Phi = lop.gaussian_dict(keys[0], m, n, normalize_atoms=True)

Let’s print the contents of the matrix

print(lop.to_matrix(Phi))
[[-0.53952553 -0.4932261   0.46508795  0.32887168  0.30488928  0.52725899
  -0.08942487 -0.95844277]
 [ 0.48259628 -0.6628074  -0.42056978 -0.27448567  0.03649064  0.80015779
   0.22029504  0.27294953]
 [ 0.06249316  0.4415227  -0.59611486 -0.46910168 -0.34679027  0.14347673
  -0.6733011   0.06830253]
 [-0.6871009   0.34996011 -0.50145921  0.77229829 -0.88625475 -0.24730527
  -0.70009922  0.047125  ]]

Check the column wise norms

print(cnb.norms_l2_cw(lop.to_matrix(Phi)))
[1. 1. 1. 1. 1. 1. 1. 1.]

Let’s print the contents of the adjoint matrix

print(lop.to_adjoint_matrix(Phi))
[[-0.53952553  0.48259628  0.06249316 -0.6871009 ]
 [-0.4932261  -0.6628074   0.4415227   0.34996011]
 [ 0.46508795 -0.42056978 -0.59611486 -0.50145921]
 [ 0.32887168 -0.27448567 -0.46910168  0.77229829]
 [ 0.30488928  0.03649064 -0.34679027 -0.88625475]
 [ 0.52725899  0.80015779  0.14347673 -0.24730527]
 [-0.08942487  0.22029504 -0.6733011  -0.70009922]
 [-0.95844277  0.27294953  0.06830253  0.047125  ]]

Let’s apply the forward operator

x = jnp.array([3, -2, 4, 8, 2, 2, 6, -1])
y = Phi.times(x)
print(y)
[  5.94539092   1.61735679 -11.34757498  -5.10351374]

Let’s apply the adjoint operator

z = Phi.trans(y)
print(z)
[  0.37032318 -10.80065612  11.40858028   2.89306098  10.32994614
   4.06291433  11.0379306   -6.27243136]

Let’s JIT compile the trans and times functions

Phi = lop.jit(Phi)

Let’s verify the forward operator after JIT compile

y = Phi.times(x)
print(y)
[  5.94539092   1.61735679 -11.34757498  -5.10351374]

Let’s verify the adjoint operator after JIT compile

z = Phi.trans(y)
print(z)
[  0.37032318 -10.80065612  11.40858028   2.89306098  10.32994614
   4.06291433  11.0379306   -6.27243136]

Column wise application on input

Number of signals

k = 5

Prepare a signal matrix

X = random.randint(keys[1], (n, k), -5, 5)
print(X)
[[-4 -5 -2 -4  2]
 [ 0  2 -2  1  3]
 [ 2 -2  4 -4 -2]
 [ 2  4 -5  1 -1]
 [ 2 -5 -2 -2 -4]
 [-4  2 -3 -1 -4]
 [-1 -1 -1  4  1]
 [-3 -4 -1 -4  1]]

Let’s apply the forward operator column wise

Y = Phi.times(X)
print(Y)
[[ 5.21151717  5.54975384  1.13780878  2.47243    -8.19423762]
 [-7.48728951 -3.88963009 -2.91612777 -2.26915609 -2.76095358]
 [-3.17949966  2.30739835 -0.17883383 -0.30940293  3.31914138]
 [ 3.06551754 12.67579867 -2.02564721  4.9074167   3.78756451]]

Let’s apply the adjoint operator column wise

Z = Phi.trans(Y)
print(Z)
[[ -8.73010144 -13.43671088  -0.64054117  -5.82025136   0.69355917]
 [  2.06116221   5.29558478   0.58377917   1.86533811   8.66256504]
 [  5.93085649  -3.51488594   2.87800132  -0.17219376  -6.52776613]
 [  7.62807664  11.59989779  -0.30588428   5.371094    -0.5688957 ]
 [ -0.29849024 -10.48404468   2.0977516   -3.57090899  -7.10687709]
 [ -4.45749663  -2.98989408  -1.25814771  -1.77009412  -6.99014907]
 [ -2.12085781 -11.7810428    0.79440614  -3.94833782  -4.76190827]
 [ -7.11129849  -5.62584597  -1.99415376  -2.77891873   7.50530159]]

Row wise application on input

Construct the operator wrapping the sensing matrix

Phi = lop.gaussian_dict(keys[0], m, n, normalize_atoms=True, axis=1)

Prepare a signal matrix with each signal along a row

XT = X.T

Let’s apply the forward operator column wise

Y = Phi.times(XT)
print(Y)
[[ 5.21151717 -7.48728951 -3.17949966  3.06551754]
 [ 5.54975384 -3.88963009  2.30739835 12.67579867]
 [ 1.13780878 -2.91612777 -0.17883383 -2.02564721]
 [ 2.47243    -2.26915609 -0.30940293  4.9074167 ]
 [-8.19423762 -2.76095358  3.31914138  3.78756451]]

Let’s apply the adjoint operator column wise

Z = Phi.trans(Y)
print(Z)
[[ -8.73010144   2.06116221   5.93085649   7.62807664  -0.29849024
   -4.45749663  -2.12085781  -7.11129849]
 [-13.43671088   5.29558478  -3.51488594  11.59989779 -10.48404468
   -2.98989408 -11.7810428   -5.62584597]
 [ -0.64054117   0.58377917   2.87800132  -0.30588428   2.0977516
   -1.25814771   0.79440614  -1.99415376]
 [ -5.82025136   1.86533811  -0.17219376   5.371094    -3.57090899
   -1.77009412  -3.94833782  -2.77891873]
 [  0.69355917   8.66256504  -6.52776613  -0.5688957   -7.10687709
   -6.99014907  -4.76190827   7.50530159]]

Square sensing matrix

Construct the operator wrapping the sensing matrix

Phi = lop.gaussian_dict(keys[0], n, normalize_atoms=False)

Let’s print the contents of the matrix

print(lop.to_matrix(Phi))
[[-5.06479786e-01 -1.76665935e-01  2.54534851e-01 -5.06364763e-01
  -1.96507134e-01  1.79333415e-01  6.42379919e-01  2.10040808e-01]
 [ 1.07240051e-01 -5.94812765e-01  1.25249038e-01  1.14459574e-01
  -2.02012833e-01 -4.09751511e-01 -4.33936884e-01 -1.34516542e-01]
 [ 1.11053768e-01  3.48672071e-01  3.68653905e-01  9.45349248e-02
  -3.29041791e-01 -1.77899332e-01 -4.04006639e-01 -1.00726911e-01]
 [ 1.33743289e-01 -5.08041866e-01 -4.99556060e-01  1.33236537e-01
  -1.68340852e-01 -2.78095509e-01 -9.48804677e-02 -3.28437869e-01]
 [-1.03963687e-03 -5.37674760e-02 -2.64978470e-01 -5.69247421e-01
  -2.17436962e-01 -3.91220044e-01 -2.19387753e-01 -2.30803981e-01]
 [ 1.71456228e-01  1.59848991e-01 -2.55230677e-01 -2.51277665e-01
   2.37702097e-02 -3.31785957e-01 -1.86295948e-01  9.79684620e-01]
 [-4.17172092e-02  2.52724150e-01  2.02487423e-01 -3.80860117e-01
  -2.98349811e-01 -1.20580149e-01  8.39526917e-01 -2.28584283e-01]
 [ 2.06869121e-01  1.04520297e+00  8.87791963e-01  2.20836838e-01
   1.17742766e-01 -1.85396223e-01  4.36565908e-03 -2.27983571e-01]]

Random matrix with orthonormal rows

Construct the operator wrapping the sensing matrix

Phi = lop.random_orthonormal_rows_dict(keys[0], m, n)

Let’s print the contents of the matrix

print(lop.to_matrix(Phi))
[[-0.48805761  0.12995918  0.4365591   0.01555415  0.05653164 -0.14781949
  -0.62155505 -0.37776645]
 [-0.46559454  0.31376094 -0.50148534  0.45550458  0.37354653  0.06455508
   0.2203263  -0.18318237]
 [ 0.34953016 -0.08652247 -0.37115078  0.21179339 -0.33253457 -0.53588747
  -0.13369294 -0.52163662]
 [-0.34327893 -0.77342544  0.15727696  0.3162613  -0.21117691  0.13210122
   0.2701341  -0.15555246]]

Check the row wise norms

print(cnb.norms_l2_rw(lop.to_matrix(Phi)))
[1. 1. 1. 1.]

Let’s apply the forward operator column wise

Y = Phi.times(X)
print(Y)
[[ 5.3156526   3.44362845  3.71437922 -0.58889084 -2.10908012]
 [ 2.58851023  4.55426428 -4.95769928  5.43997158 -1.15770268]
 [ 1.46024812  2.47991672 -0.14152435  2.96448472  3.78835967]
 [ 1.56595683  2.79219775  1.19267863  2.27984255 -3.206765  ]]

Let’s apply the adjoint operator column wise

Z = Phi.trans(Y)
print(Z)
[[-3.82670014 -3.8928266   0.03655821 -1.9918535   3.99333346]
 [ 0.16549613 -0.49764393 -1.98301766 -0.38946395  1.51507918]
 [ 0.72681332 -1.26182404  4.3478673  -3.72685574 -2.2505701 ]
 [ 2.06628118  3.53634503 -1.85325656  3.81765653 -0.77196994]
 [ 0.45115468  0.48159786 -1.84675587  0.53154779 -1.13425141]
 [-1.19431947 -1.17513811 -0.63570689 -0.84923291 -2.21672446]
 [-2.50586037 -0.71426    -3.05989872  1.7841294  -0.31689645]
 [-3.48755201 -3.86309685 -0.60670455 -2.67506262 -0.46851652]]

The frame operator

F = lop.frame(Phi)

Check that it’s close to identity matrix

print(jnp.round(lop.to_matrix(F),2))
[[ 1. -0.  0.  0.]
 [-0.  1.  0.  0.]
 [ 0.  0.  1. -0.]
 [ 0.  0. -0.  1.]]

Random orthonormal basis operator

Construct the operator wrapping the sensing matrix

Phi = lop.random_onb_dict(keys[0], n)

Let’s print the contents of the matrix

print(lop.to_matrix(Phi))
[[-0.83021008  0.10099603  0.1924291  -0.18709837 -0.05440769 -0.42034524
   0.05004106  0.21531801]
 [ 0.17578544 -0.49159467  0.73866229 -0.34880127  0.22399819  0.03967139
  -0.03076634  0.08633818]
 [ 0.1820368   0.2083679   0.2342978  -0.08887811 -0.75177769 -0.00617792
  -0.49250229  0.22996581]
 [ 0.21922894 -0.43921983 -0.12485129  0.04649406 -0.44739928 -0.55798428
   0.46709786 -0.10762012]
 [-0.00170415 -0.03951229 -0.294542   -0.6724019   0.03554208 -0.17193851
  -0.34738298 -0.55502752]
 [ 0.28104713  0.04021688 -0.36937749 -0.38856286  0.23274741 -0.19465774
   0.00554408  0.73476564]
 [-0.06838189  0.20710371  0.01707164 -0.47829639 -0.29978828  0.51519234
   0.60668761 -0.00266108]
 [ 0.33909513  0.6824461   0.35125433 -0.0373899   0.19031866 -0.4212764
   0.21687527 -0.18377793]]

Check the row wise norms

print(cnb.norms_l2_rw(lop.to_matrix(Phi)))
[1. 1. 1. 1. 1. 1. 1. 1.]

Check the column wise norms

print(cnb.norms_l2_cw(lop.to_matrix(Phi)))
[1. 1. 1. 1. 1. 1. 1. 1.]

Let’s apply the forward operator column wise

Y = Phi.times(X)
print(Y)
[[ 4.20807229  1.73982564  4.26812837  2.33307439  0.60917887]
 [ 0.13764291 -6.08988076  4.70769159 -5.45427267 -3.25084291]
 [-2.11354663  2.00161547  2.38539797 -2.92598768  3.37874579]
 [ 0.15927087 -0.45449504  1.91737998  2.9814184   3.7050189 ]
 [ 0.84423256 -0.12512179  3.61341624  1.4045029   0.78271592]
 [-3.60578444 -6.63795754 -0.79905496 -3.18274785  2.3980138 ]
 [-3.90797221  0.74207081  0.63229793  2.4558273   0.67111076]
 [ 1.67154977 -2.45856137  0.39897891 -0.4730898   3.0173381 ]]

Let’s apply the adjoint operator column wise

Z = Phi.trans(Y)
print(Z)
[[-4. -5. -2. -4.  2.]
 [ 0.  2. -2.  1.  3.]
 [ 2. -2.  4. -4. -2.]
 [ 2.  4. -5.  1. -1.]
 [ 2. -5. -2. -2. -4.]
 [-4.  2. -3. -1. -4.]
 [-1. -1. -1.  4.  1.]
 [-3. -4. -1. -4.  1.]]

The gram operator

G = lop.gram(Phi)

Check that it’s close to identity matrix

print(jnp.round(lop.to_matrix(G),2))
[[ 1.  0.  0.  0.  0.  0.  0. -0.]
 [ 0.  1.  0. -0.  0.  0. -0. -0.]
 [ 0.  0.  1. -0. -0.  0.  0.  0.]
 [ 0. -0. -0.  1.  0. -0. -0.  0.]
 [ 0.  0. -0.  0.  1. -0. -0. -0.]
 [ 0.  0.  0. -0. -0.  1. -0.  0.]
 [ 0. -0.  0. -0. -0. -0.  1.  0.]
 [-0. -0.  0.  0. -0.  0.  0.  1.]]

The frame operator

F = lop.frame(Phi)

Check that it’s close to identity matrix

print(jnp.round(lop.to_matrix(F),2))
[[ 1. -0. -0.  0.  0.  0.  0.  0.]
 [-0.  1.  0.  0.  0. -0.  0. -0.]
 [-0.  0.  1. -0.  0.  0. -0. -0.]
 [ 0.  0. -0.  1.  0.  0. -0.  0.]
 [ 0.  0.  0.  0.  1. -0.  0.  0.]
 [ 0. -0.  0.  0. -0.  1.  0.  0.]
 [ 0.  0. -0. -0.  0.  0.  1.  0.]
 [ 0. -0. -0.  0.  0.  0.  0.  1.]]

Total running time of the script: (0 minutes 2.337 seconds)

Gallery generated by Sphinx-Gallery