Skip to content

Commit a92eab8

Browse files
committed
Added BPDN tutorial
1 parent dfcbb89 commit a92eab8

File tree

1 file changed

+35
-2
lines changed

1 file changed

+35
-2
lines changed

tutorials/basispursuit.py

Lines changed: 35 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,11 @@
88
\argmin_\mathbf{x} \|\mathbf{x}\|_1 \; \text{s.t.} \; \mathbf{Ax} = \mathbf{y}
99
1010
where the operator :math:`\mathbf{A}` is of size :math:`N \times M`, and generally :math:`N<M`.
11-
Note that this problem is similar to more general L1-regularized inversion, but it presents a stricter condition on the
12-
data term which must be satisfied exactly.
11+
Note that this problem is similar to more general L1-regularized inversion, but it presents a stricter condition
12+
on the data term which must be satisfied exactly. Similarly, we can also consider the Basis Pursuit Denoise problem
13+
14+
.. math::
15+
\argmin_\mathbf{x} \|\mathbf{x}\|_1 \; \text{s.t.} \; \|\mathbf{Ax} - \mathbf{y}\|_2 < \epsilon
1316
1417
"""
1518
import numpy as np
@@ -65,3 +68,33 @@
6568
# We can observe how even after few iterations, despite the solution is not
6669
# yet converged the data reconstruction is perfect. This is consequence of the
6770
# fact that for the Basis Pursuit problem the data term is a hard constraint.
71+
# Finally let's consider the case with a soft constraint.
72+
73+
f = pyproximal.L1()
74+
g = pyproximal.proximal.EuclideanBall(y, .1)
75+
76+
L = np.real((Aop.H @ Aop).eigs(1))[0]
77+
tau = .99
78+
mu = tau / L
79+
80+
xinv_early = pyproximal.optimization.primaldual.PrimalDual(f, g, Aop, np.zeros_like(x),
81+
tau, mu, niter=10)
82+
83+
xinv = pyproximal.optimization.primaldual.PrimalDual(f, g, Aop, np.zeros_like(x),
84+
tau, mu, niter=1000)
85+
86+
fig, axs = plt.subplots(1, 2, figsize=(12, 3))
87+
axs[0].plot(x, 'k')
88+
axs[0].plot(xinv_early, '--r')
89+
axs[0].plot(xinv, '--b')
90+
axs[0].set_title('Model')
91+
axs[1].plot(y, 'k', label='True')
92+
axs[1].plot(Aop * xinv_early, '--r', label='Early Inv')
93+
axs[1].plot(Aop * xinv, '--b', label='Inv (Res=%.2f)' % np.linalg.norm(y - Aop @ xinv))
94+
axs[1].set_title('Data')
95+
axs[1].legend()
96+
plt.tight_layout()
97+
98+
###############################################################################
99+
# Note that at convergence the norm of the residual of the solution adheres
100+
# to the EuclideanBall constraint

0 commit comments

Comments
 (0)