Merge pull request #95 from mrava87/main

Added BPDN tutorial

Author: mrava87

tutorials/basispursuit.py

@@ -8,8 +8,11 @@
     \argmin_\mathbf{x} \|\mathbf{x}\|_1 \; \text{s.t.} \; \mathbf{Ax} = \mathbf{y}
 
 where the operator :math:`\mathbf{A}` is of size :math:`N \times M`, and generally :math:`N<M`.
-Note that this problem is similar to more general L1-regularized inversion, but it presents a stricter condition on the
-data term which must be satisfied exactly.
+Note that this problem is similar to more general L1-regularized inversion, but it presents a stricter condition
+on the data term which must be satisfied exactly. Similarly, we can also consider the Basis Pursuit Denoise problem
+
+.. math::
+    \argmin_\mathbf{x} \|\mathbf{x}\|_1 \; \text{s.t.} \;  \|\mathbf{Ax} - \mathbf{y}\|_2 < \epsilon
 
 """
 import numpy as np
@@ -65,3 +68,33 @@
 # We can observe how even after few iterations, despite the solution is not
 # yet converged the data reconstruction is perfect. This is consequence of the
 # fact that for the Basis Pursuit problem the data term is a hard constraint.
+# Finally let's consider the case with a soft constraint.
+
+f = pyproximal.L1()
+g = pyproximal.proximal.EuclideanBall(y, .1)
+
+L = np.real((Aop.H @ Aop).eigs(1))[0]
+tau = .99
+mu = tau / L
+
+xinv_early = pyproximal.optimization.primaldual.PrimalDual(f, g, Aop, np.zeros_like(x),
+                                                           tau, mu, niter=10)
+
+xinv = pyproximal.optimization.primaldual.PrimalDual(f, g, Aop, np.zeros_like(x),
+                                                     tau, mu, niter=1000)
+
+fig, axs = plt.subplots(1, 2, figsize=(12, 3))
+axs[0].plot(x, 'k')
+axs[0].plot(xinv_early, '--r')
+axs[0].plot(xinv, '--b')
+axs[0].set_title('Model')
+axs[1].plot(y, 'k', label='True')
+axs[1].plot(Aop * xinv_early, '--r', label='Early Inv')
+axs[1].plot(Aop * xinv, '--b', label='Inv (Res=%.2f)' % np.linalg.norm(y - Aop @ xinv))
+axs[1].set_title('Data')
+axs[1].legend()
+plt.tight_layout()
+
+###############################################################################
+# Note that at convergence the norm of the residual of the solution adheres
+# to the EuclideanBall constraint