Note
Click here 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.sparse as crs
# 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))
Out:
[[-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))
Out:
[[-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)
Out:
[-24. -2. -18. 16.]
Let’s apply the adjoint operator
z = Phi.trans(y)
print(z)
Out:
[ 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)
Out:
[-24. -2. -18. 16.]
Let’s verify the adjoint operator after JIT compile
z = Phi.trans(y)
print(z)
Out:
[ 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))
Out:
[[-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(crs.norms_l2_cw(lop.to_matrix(Phi)))
Out:
[1. 1. 1. 1. 1. 1. 1. 1.]
Let’s print the contents of the adjoint matrix
print(lop.to_adjoint_matrix(Phi))
Out:
[[-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)
Out:
[-12. -1. -9. 8.]
Let’s apply the adjoint operator
z = Phi.trans(y)
print(z)
Out:
[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)
Out:
[-12. -1. -9. 8.]
Let’s verify the adjoint operator after JIT compile
z = Phi.trans(y)
print(z)
Out:
[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)
Out:
[[-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)
Out:
[[ 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)
Out:
[[ -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)
Out:
[[ 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)
Out:
[[ -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))
Out:
[[ 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))
Out:
[[-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))
Out:
[[-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)
Out:
[ 5.79912122 1.61042927 -12.54007687 -3.16525231]
Let’s apply the adjoint operator
z = Phi.trans(y)
print(z)
Out:
[ -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)
Out:
[ 5.79912122 1.61042927 -12.54007687 -3.16525231]
Let’s verify the adjoint operator after JIT compile
z = Phi.trans(y)
print(z)
Out:
[ -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))
Out:
[[-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(crs.norms_l2_cw(lop.to_matrix(Phi)))
Out:
[1. 1. 1. 1. 1. 1. 1. 1.]
Let’s print the contents of the adjoint matrix
print(lop.to_adjoint_matrix(Phi))
Out:
[[-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)
Out:
[ 5.94539092 1.61735679 -11.34757498 -5.10351374]
Let’s apply the adjoint operator
z = Phi.trans(y)
print(z)
Out:
[ 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)
Out:
[ 5.94539092 1.61735679 -11.34757498 -5.10351374]
Let’s verify the adjoint operator after JIT compile
z = Phi.trans(y)
print(z)
Out:
[ 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)
Out:
[[-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)
Out:
[[ 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)
Out:
[[ -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)
Out:
[[ 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)
Out:
[[ -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))
Out:
[[-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))
Out:
[[-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(crs.norms_l2_rw(lop.to_matrix(Phi)))
Out:
[1. 1. 1. 1.]
Let’s apply the forward operator column wise
Y = Phi.times(X)
print(Y)
Out:
[[ 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)
Out:
[[-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))
Out:
[[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))
Out:
[[-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(crs.norms_l2_rw(lop.to_matrix(Phi)))
Out:
[1. 1. 1. 1. 1. 1. 1. 1.]
Check the column wise norms
print(crs.norms_l2_cw(lop.to_matrix(Phi)))
Out:
[1. 1. 1. 1. 1. 1. 1. 1.]
Let’s apply the forward operator column wise
Y = Phi.times(X)
print(Y)
Out:
[[ 5.57873951 5.82895517 1.38843001 5.56891952 0.47398041]
[-0.96530474 0.28435132 -4.57847575 2.08946428 -1.10494 ]
[-2.30043511 4.50455341 -0.90283847 4.90589573 4.66325372]
[ 0.76117474 -1.25261796 4.10008637 -2.70812388 2.57964641]
[ 0.68707703 -1.05906074 -2.55908778 -1.2681856 0.06456515]
[ 3.84565902 0.88224872 -0.24904421 0.47313045 -1.33724135]
[-0.86566444 5.82881512 -3.07918607 -1.32579627 -1.42524853]
[-0.25282498 -1.79089772 -2.71876006 -0.79350088 4.28138095]]
Let’s apply the adjoint operator column wise
Z = Phi.trans(Y)
print(Z)
Out:
[[-4.00000000e+00 -5.00000000e+00 -2.00000000e+00 -4.00000000e+00
2.00000000e+00]
[-2.77555756e-17 2.00000000e+00 -2.00000000e+00 1.00000000e+00
3.00000000e+00]
[ 2.00000000e+00 -2.00000000e+00 4.00000000e+00 -4.00000000e+00
-2.00000000e+00]
[ 2.00000000e+00 4.00000000e+00 -5.00000000e+00 1.00000000e+00
-1.00000000e+00]
[ 2.00000000e+00 -5.00000000e+00 -2.00000000e+00 -2.00000000e+00
-4.00000000e+00]
[-4.00000000e+00 2.00000000e+00 -3.00000000e+00 -1.00000000e+00
-4.00000000e+00]
[-1.00000000e+00 -1.00000000e+00 -1.00000000e+00 4.00000000e+00
1.00000000e+00]
[-3.00000000e+00 -4.00000000e+00 -1.00000000e+00 -4.00000000e+00
1.00000000e+00]]
The gram operator
G = lop.gram(Phi)
Check that it’s close to identity matrix
print(jnp.round(lop.to_matrix(G),2))
Out:
[[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))
Out:
[[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 3.077 seconds)
