Note
Go to the end to download the full example code
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
Unnormalized atoms / columns¶
Construct the operator wrapping the sensing matrix
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
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
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
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
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
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
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.243 seconds)
