Compressive Sensing with factorized operator¶
This notebooks consider a compressive sensing problem:
$$ \arg min_\mathbf{x} = \frac{1}{2} \|\mathbf{Ax}-\mathbf{y}\|_2^2 + \lambda \|\mathbf{x}\|_1 $$where the matrix $\mathbf{A}$ is a randomized sampling matrix.
We investigate here the scenario where the matrix $A$ is dense and can be explicitely factorized via Cholesky factorization and used in subsequent operations when computing the proximal of the L2 operator. We will see that this drammatically reduce the cost of the L2 operator, and the overall solution of the problem at hand.
In [1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pylops
import pyproximal
In [2]:
n, m = 40, 100
mava = 10
# model
x = np.zeros(m)
iava = np.random.permutation(np.arange(m))[:mava]
x[iava] = np.random.normal(0, 1, mava)
# operator
A = np.random.normal(0, 1, (n, m))
Aop = pylops.MatrixMult(A)
y = Aop * x
In [3]:
fig, axs = plt.subplots(1, 2, figsize=(12, 3))
axs[0].plot(x, 'k')
axs[0].set_title('x')
axs[1].plot(y, 'k')
axs[1].set_title('y');
In [4]:
f = pyproximal.L2(Aop, y, niter=20)
f1 = pyproximal.L2(Aop, y, densesolver='scipy')
f2 = pyproximal.L2(Aop, y, densesolver='factorize')
g = pyproximal.L1()
ATA = Aop.H * Aop
ATA.explicit = False
L = np.abs(ATA.eigs(1)[0])
niter = 1000
tau= 0.95/L
xinv = pyproximal.optimization.primal.ADMM(f, g, np.zeros_like(x),
tau, niter=niter, show=True)[0]
xinv1 = pyproximal.optimization.primal.ADMM(f1, g, np.zeros_like(x),
tau, niter=niter, show=True)[0]
xinv2 = pyproximal.optimization.primal.ADMM(f2, g, np.zeros_like(x),
tau, niter=niter, show=True)[0]
fig, axs = plt.subplots(2, 1, figsize=(16, 13))
axs[0].plot(x, 'k')
axs[0].plot(xinv, '--r')
axs[0].plot(xinv1, '--b')
axs[0].plot(xinv2, '--g')
axs[0].set_title('x')
axs[1].plot(y, 'k')
axs[1].plot(Aop @ xinv, '--r')
axs[1].plot(Aop @ xinv1, '--b')
axs[1].plot(Aop @ xinv2, '--g')
axs[1].set_title('y');
ADMM
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.L1.L1'>
tau = 3.584829e-03 niter = 1000
Itn x[0] f g J = f + g
1 2.70646e-02 1.121e+02 4.317e+00 1.164e+02
2 2.92837e-02 5.941e+01 6.658e+00 6.607e+01
3 2.58756e-02 3.308e+01 8.451e+00 4.153e+01
4 1.81746e-02 2.002e+01 9.655e+00 2.968e+01
5 9.08583e-03 1.297e+01 1.050e+01 2.348e+01
6 -2.08513e-04 8.912e+00 1.109e+01 2.000e+01
7 -8.68047e-03 6.421e+00 1.150e+01 1.792e+01
8 -4.00978e-03 4.823e+00 1.178e+01 1.660e+01
9 -5.81808e-03 3.727e+00 1.197e+01 1.570e+01
10 -7.15890e-03 2.961e+00 1.211e+01 1.507e+01
101 2.02745e-05 4.998e-01 1.026e+01 1.076e+01
201 -1.23418e-05 4.613e-01 9.332e+00 9.793e+00
301 1.02329e-05 4.280e-01 8.857e+00 9.284e+00
401 6.96832e-06 4.216e-01 8.587e+00 9.009e+00
501 -1.51011e-06 4.001e-01 8.403e+00 8.803e+00
601 -5.56403e-06 3.538e-01 8.302e+00 8.656e+00
701 1.24165e-06 2.892e-01 8.273e+00 8.562e+00
801 3.41502e-06 2.364e-01 8.281e+00 8.517e+00
901 5.16943e-07 2.023e-01 8.299e+00 8.501e+00
992 3.28505e-07 1.858e-01 8.313e+00 8.499e+00
993 3.24217e-07 1.858e-01 8.314e+00 8.499e+00
994 3.19981e-07 1.857e-01 8.314e+00 8.499e+00
995 3.15795e-07 1.856e-01 8.314e+00 8.499e+00
996 3.11660e-07 1.855e-01 8.314e+00 8.499e+00
997 3.07575e-07 1.854e-01 8.314e+00 8.499e+00
998 3.03539e-07 1.853e-01 8.314e+00 8.499e+00
999 2.99553e-07 1.852e-01 8.314e+00 8.499e+00
1000 2.95616e-07 1.851e-01 8.314e+00 8.499e+00
Total time (s) = 0.20
---------------------------------------------------------
ADMM
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.L1.L1'>
tau = 3.584829e-03 niter = 1000
Itn x[0] f g J = f + g
1 2.70646e-02 1.121e+02 4.317e+00 1.164e+02
2 2.92837e-02 5.941e+01 6.658e+00 6.607e+01
3 2.58756e-02 3.308e+01 8.451e+00 4.153e+01
4 1.81746e-02 2.002e+01 9.655e+00 2.968e+01
5 9.08583e-03 1.297e+01 1.050e+01 2.348e+01
6 -2.08513e-04 8.912e+00 1.109e+01 2.000e+01
7 -8.68047e-03 6.421e+00 1.150e+01 1.792e+01
8 -4.00978e-03 4.823e+00 1.178e+01 1.660e+01
9 -5.81808e-03 3.727e+00 1.197e+01 1.570e+01
10 -7.15890e-03 2.961e+00 1.211e+01 1.507e+01
101 2.02745e-05 4.998e-01 1.026e+01 1.076e+01
201 -1.23418e-05 4.613e-01 9.332e+00 9.793e+00
301 1.02329e-05 4.280e-01 8.857e+00 9.284e+00
401 6.96832e-06 4.216e-01 8.587e+00 9.009e+00
501 -1.51011e-06 4.001e-01 8.403e+00 8.803e+00
601 -5.56403e-06 3.538e-01 8.302e+00 8.656e+00
701 1.24165e-06 2.892e-01 8.273e+00 8.562e+00
801 3.41502e-06 2.364e-01 8.281e+00 8.517e+00
901 5.16943e-07 2.023e-01 8.299e+00 8.501e+00
992 3.28505e-07 1.858e-01 8.313e+00 8.499e+00
993 3.24217e-07 1.858e-01 8.314e+00 8.499e+00
994 3.19981e-07 1.857e-01 8.314e+00 8.499e+00
995 3.15795e-07 1.856e-01 8.314e+00 8.499e+00
996 3.11660e-07 1.855e-01 8.314e+00 8.499e+00
997 3.07575e-07 1.854e-01 8.314e+00 8.499e+00
998 3.03539e-07 1.853e-01 8.314e+00 8.499e+00
999 2.99553e-07 1.852e-01 8.314e+00 8.499e+00
1000 2.95616e-07 1.851e-01 8.314e+00 8.499e+00
Total time (s) = 0.18
---------------------------------------------------------
ADMM
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.L1.L1'>
tau = 3.584829e-03 niter = 1000
Itn x[0] f g J = f + g
1 2.70646e-02 1.121e+02 4.317e+00 1.164e+02
2 2.92837e-02 5.941e+01 6.658e+00 6.607e+01
3 2.58756e-02 3.308e+01 8.451e+00 4.153e+01
4 1.81746e-02 2.002e+01 9.655e+00 2.968e+01
5 9.08583e-03 1.297e+01 1.050e+01 2.348e+01
6 -2.08513e-04 8.912e+00 1.109e+01 2.000e+01
7 -8.68047e-03 6.421e+00 1.150e+01 1.792e+01
8 -4.00978e-03 4.823e+00 1.178e+01 1.660e+01
9 -5.81808e-03 3.727e+00 1.197e+01 1.570e+01
10 -7.15890e-03 2.961e+00 1.211e+01 1.507e+01
101 2.02745e-05 4.998e-01 1.026e+01 1.076e+01
201 -1.23418e-05 4.613e-01 9.332e+00 9.793e+00
301 1.02329e-05 4.280e-01 8.857e+00 9.284e+00
401 6.96832e-06 4.216e-01 8.587e+00 9.009e+00
501 -1.51011e-06 4.001e-01 8.403e+00 8.803e+00
601 -5.56403e-06 3.538e-01 8.302e+00 8.656e+00
701 1.24165e-06 2.892e-01 8.273e+00 8.562e+00
801 3.41502e-06 2.364e-01 8.281e+00 8.517e+00
901 5.16943e-07 2.023e-01 8.299e+00 8.501e+00
992 3.28505e-07 1.858e-01 8.313e+00 8.499e+00
993 3.24217e-07 1.858e-01 8.314e+00 8.499e+00
994 3.19981e-07 1.857e-01 8.314e+00 8.499e+00
995 3.15795e-07 1.856e-01 8.314e+00 8.499e+00
996 3.11660e-07 1.855e-01 8.314e+00 8.499e+00
997 3.07575e-07 1.854e-01 8.314e+00 8.499e+00
998 3.03539e-07 1.853e-01 8.314e+00 8.499e+00
999 2.99553e-07 1.852e-01 8.314e+00 8.499e+00
1000 2.95616e-07 1.851e-01 8.314e+00 8.499e+00
Total time (s) = 0.05
---------------------------------------------------------
In [5]:
%timeit -n 5 -r 3 pyproximal.optimization.primal.ADMM(f, g, np.zeros_like(x), tau, niter=niter)
%timeit -n 5 -r 3 pyproximal.optimization.primal.ADMM(f1, g, np.zeros_like(x), tau, niter=niter)
%timeit -n 5 -r 3 pyproximal.optimization.primal.ADMM(f2, g, np.zeros_like(x), tau, niter=niter)
172 ms ± 4.34 ms per loop (mean ± std. dev. of 3 runs, 5 loops each) 163 ms ± 1.18 ms per loop (mean ± std. dev. of 3 runs, 5 loops each) 33.6 ms ± 240 µs per loop (mean ± std. dev. of 3 runs, 5 loops each)
In [6]:
np.linalg.norm(x-xinv), np.linalg.norm(x-xinv1), np.linalg.norm(x-xinv2)
Out[6]:
(0.2038399703562013, 0.2038399703562013, 0.203839970356175)