Compressed Sensing & Proximal Methods
=====================================

Introduction
-------------

Compressed sensing is a revolutionary signal processing technique that exploits the sparsity of signals to reconstruct them from far fewer measurements than those required by the classical Nyquist sampling theorem :cite:`donoho_2006_cs,candes_romberg_tao_2006`. In radio interferometric imaging, this is particularly powerful because astronomical sources are often sparse in appropriate bases :cite:`wiaux_2009_cs_ri,carrillo_2014_purify`.

This document provides a comprehensive introduction to compressed sensing theory and its implementation through proximal optimization methods. We begin with the mathematical foundations, explain why traditional optimization approaches fail, and then introduce proximal operators as the solution. Finally, we present practical algorithms (FISTA and SDMM) for solving compressed sensing problems.

Mathematical Foundation of Compressed Sensing
---------------------------------------------

**The Nyquist Sampling Limit:**

The traditional Nyquist sampling theorem states that to perfectly reconstruct a signal, we must sample it at least twice the highest frequency present in the signal. For a signal with bandwidth :math:`B`, this requires at least :math:`2B` samples per second. In radio interferometry, this translates to needing complete coverage of the :math:`(u,v)` plane to reconstruct the image without aliasing.

**Breaking the Nyquist Barrier:**

Compressed sensing challenges this fundamental limit by showing that if a signal is sparse in some basis, we can reconstruct it from far fewer measurements than Nyquist requires :cite:`donoho_2006_cs,candes_romberg_tao_2006`. Instead of needing :math:`n` measurements for an :math:`n`-dimensional signal, we only need roughly :math:`k \log(n/k)` measurements, where :math:`k` is the number of non-zero coefficients in the sparse representation.

This dramatic reduction in measurement requirements has profound implications for radio interferometry, where complete :math:`(u,v)` coverage is often impossible due to practical constraints.

**The Sparsity Assumption:**

The key insight of compressed sensing is that many signals can be represented as a linear combination of basis functions with sparse coefficients:

.. math::

   \mathbf{s} = \mathbf{\Phi} \mathbf{x}

where:

- :math:`\mathbf{s}` is the original signal (image)
- :math:`\mathbf{\Phi}` is the basis matrix (e.g., Fourier, wavelet, etc.)
- :math:`\mathbf{x}` is the sparse coefficient vector

If the image is already sparse in the pixel domain, then :math:`\mathbf{\Phi}` is simply the identity matrix (unitary).

**Formulating the Compressed Sensing Problem:**

In traditional linear inverse problems, we solve:

.. math::

   \mathbf{A}\mathbf{s} + \mathbf{\epsilon} = \mathbf{b}

where :math:`\mathbf{A}` is the measurement matrix and :math:`\mathbf{b}` are the measurements.

For compressed sensing, this becomes:

.. math::

   \mathbf{A}\mathbf{\Phi}\mathbf{x} + \mathbf{\epsilon} = \mathbf{y}

where:

- :math:`\mathbf{A}` is the measurement matrix (e.g., sampling pattern)
- :math:`\mathbf{\Phi}` is the sparsity basis
- :math:`\mathbf{x}` is the sparse coefficient vector
- :math:`\mathbf{y}` are the compressed measurements

This formulation allows us to work directly with the sparse coefficients :math:`\mathbf{x}` rather than the full signal :math:`\mathbf{s}`, enabling reconstruction from fewer measurements.

From Theory to Optimization
----------------------------

**The Sparsity Optimization Challenge:**

To reconstruct the sparse signal, we need to solve an optimization problem. In principle, we would like to minimize the L0 norm (the number of non-zero elements):

.. math::

   \|\mathbf{x}\|_0 = \text{number of non-zero elements in } \mathbf{x}

However, this leads to a computationally intractable NP-hard combinatorial problem that becomes prohibitive for large signals.

**The L1 Relaxation Solution:**

The L1 norm provides a convex relaxation of the L0 norm:

.. math::

   \|\mathbf{x}\|_1 = \sum_i |x_i|

This relaxation has several key advantages :cite:`donoho_2006_cs`:

- **Convex**: Makes the optimization problem tractable
- **Sparsity-promoting**: Under certain conditions, L1 minimization yields sparse solutions
- **Computationally efficient**: Can be solved using convex optimization methods

**The Complete Optimization Problem:**

We solve the following convex optimization problem:

.. math::

   \min_{\mathbf{x}} \|\mathbf{A}\mathbf{\Phi}\mathbf{x} - \mathbf{y}\|_2^2 + \lambda \|\mathbf{x}\|_1

where:

- The first term is the data fidelity (least squares)
- The second term is the L1 sparsity penalty
- :math:`\lambda` controls the trade-off between fidelity and sparsity

**The Differentiability Challenge:**

While the L1 norm is convex and promotes sparsity, it is not differentiable at zero, which poses challenges for gradient-based optimization methods. This is where proximal operators become essential.

Theoretical Guarantees
----------------------

**When Does Compressed Sensing Work?**

For compressed sensing to succeed, three key conditions must be satisfied:

**1. Incoherence Condition:**
The measurement matrix :math:`\mathbf{A}` must be incoherent with respect to the sparsity basis :math:`\mathbf{\Phi}`. Incoherence means that the measurement vectors cannot be too aligned with the sparsity basis vectors.

Mathematically, the coherence :math:`\mu` between :math:`\mathbf{A}` and :math:`\mathbf{\Phi}` is defined as:

.. math::

   \mu(\mathbf{A}, \mathbf{\Phi}) = \max_{i,j} |\langle \mathbf{a}_i, \mathbf{\phi}_j \rangle|

where :math:`\mathbf{a}_i` are the rows of :math:`\mathbf{A}` and :math:`\mathbf{\phi}_j` are the columns of :math:`\mathbf{\Phi}`.

For compressed sensing to work, we need :math:`\mu(\mathbf{A}, \mathbf{\Phi})` to be small (typically :math:`\mu \ll 1`).

**2. Sparsity Level Requirement:**
The number of measurements :math:`p` must satisfy:

.. math::

   p \sim O(k \log(n/k))

where:

- :math:`k` is the number of non-zero elements in :math:`\mathbf{x}`
- :math:`n` is the dimension of :math:`\mathbf{x}`

This means we need roughly :math:`k \log(n/k)` measurements to recover a :math:`k`-sparse signal.

**3. Restricted Isometry Property (RIP):**
The matrix :math:`\mathbf{A}\mathbf{\Phi}` must satisfy the RIP, which means it acts approximately like a unitary matrix when restricted to sparse vectors.

For a :math:`k`-sparse vector :math:`\mathbf{x}`, the RIP constant :math:`\delta_k` satisfies:

.. math::

   (1 - \delta_k)\|\mathbf{x}\|_2^2 \leq \|\mathbf{A}\mathbf{\Phi}\mathbf{x}\|_2^2 \leq (1 + \delta_k)\|\mathbf{x}\|_2^2

For successful recovery, we typically need :math:`\delta_{2k} < 0.414`.

**Application to Radio Interferometry:**

In radio interferometry, these conditions translate to:

  - :math:`\mathbf{A}` represents the sampling pattern in the :math:`(u,v)` plane
  - :math:`\mathbf{\Phi}` could be the Fourier basis (for point sources) or wavelet basis (for extended sources)
  - :math:`\mathbf{y}` are the measured visibilities
  - The sparsity assumption is often valid for astronomical sources

The incoherence condition is particularly important: random sampling patterns in the :math:`(u,v)` plane are naturally incoherent with Fourier bases, making compressed sensing well-suited for radio interferometry.

Proximal Operators: The Key to Practical Implementation
-------------------------------------------------------

As we saw above, the L1 norm is not differentiable at zero, which prevents the use of standard gradient-based optimization methods. Proximal operators provide the solution to this challenge, making compressed sensing optimization tractable.

**Definition:**

The proximal operator of a function :math:`g` at point :math:`\mathbf{x}` is defined as:

.. math::

   \text{prox}_{\lambda g}(\mathbf{x}) = \arg\min_{\mathbf{z}} \left\{ \lambda g(\mathbf{z}) + \frac{1}{2}\|\mathbf{z} - \mathbf{x}\|_2^2 \right\}

where :math:`\lambda > 0` is a parameter controlling the strength of the penalty.

**Key Properties:**

- **Nonexpansive**: :math:`\|\text{prox}_{\lambda g}(\mathbf{x}) - \text{prox}_{\lambda g}(\mathbf{y})\|_2 \leq \|\mathbf{x} - \mathbf{y}\|_2`
- **Firmly nonexpansive**: Provides strong convergence guarantees
- **Computable**: For many functions, the proximal operator has a closed-form solution

**Common Proximal Operators:**

**1. L1 Norm (Soft Thresholding):**
For :math:`g(\mathbf{x}) = \|\mathbf{x}\|_1`, the proximal operator is:

.. math::

   \text{prox}_{\lambda \|\cdot\|_1}(\mathbf{x})_i = \text{sign}(x_i) \cdot \max(|x_i| - \lambda, 0)

This is known as **soft thresholding** and promotes sparsity by shrinking small coefficients toward zero.

**2. L2 Norm Squared:**
For :math:`g(\mathbf{x}) = \frac{1}{2}\|\mathbf{x}\|_2^2`, the proximal operator is:

.. math::

   \text{prox}_{\lambda \|\cdot\|_2^2}(\mathbf{x}) = \frac{\mathbf{x}}{1 + \lambda}

**3. Indicator Function:**
For a convex set :math:`C`, the indicator function :math:`\iota_C(\mathbf{x}) = 0` if :math:`\mathbf{x} \in C` and :math:`+\infty` otherwise. The proximal operator is:

.. math::

   \text{prox}_{\iota_C}(\mathbf{x}) = \arg\min_{\mathbf{z} \in C} \|\mathbf{z} - \mathbf{x}\|_2^2

This is the **projection** onto the set :math:`C`.

**Why Proximal Operators Work:**

Proximal operators provide a way to "softly" enforce constraints. Instead of hard constraints (which can be difficult to handle), they penalize violations while allowing some flexibility. This makes optimization problems more tractable while still achieving the desired sparsity or other properties.

**The Connection to Compressed Sensing:**

For the L1 norm, the proximal operator implements soft-thresholding, which naturally promotes sparsity by shrinking small coefficients toward zero while preserving larger ones. This elegant solution allows us to handle the non-differentiability of the L1 norm while maintaining the sparsity-promoting properties essential for compressed sensing reconstruction.

Proximal Optimization Algorithms
---------------------------------

Now that we understand the mathematical foundation and the role of proximal operators, we can implement efficient optimization algorithms. Pyralysis provides two main proximal optimization methods:

+---------------------+-----------------------------------------------+
| Optimizer           | Description                                   |
+=====================+===============================================+
| FISTA               | Fast Iterative Shrinkage-Thresholding, O(1/k²)|
+---------------------+-----------------------------------------------+
| MFISTA              | Monotone FISTA, ensures objective decrease    |
+---------------------+-----------------------------------------------+
| SDMM                | Simultaneous Direction Method of Multipliers  |
+---------------------+-----------------------------------------------+

FISTA (Fast Iterative Shrinkage-Thresholding Algorithm)
-------------------------------------------------------

FISTA is a proximal gradient method that achieves O(1/k²) convergence rate compared to O(1/k) for ISTA, where :math:`k` is the number of iterations :cite:`beck_teboulle_2009_fista`. It's particularly effective for problems with exactly one non-differentiable term.

**Mathematical Foundation:**

FISTA solves problems of the form:

.. math::

   \arg\min_x f(x) + \lambda g(x)

where:

- :math:`f(x)` is smooth (data fidelity term)
- :math:`g(x)` is non-differentiable (regularization term)
- :math:`\lambda` is the regularization parameter

**Algorithm Steps:**

1. Compute gradient of smooth term: :math:`\nabla f(x)`
2. Apply line search for step size
3. Apply proximal operator: :math:`\text{prox}_{\lambda g}(x)`
4. Update acceleration parameters

**Proximal Operator for L1 Regularization:**
For L1 sparsity, the proximal operator implements soft-thresholding:

.. math::

   \text{prox}_{\lambda \|\cdot\|_1}(x) = \text{sign}(x) \cdot \max(|x| - \lambda, 0)

**FISTA Variants:**

**Standard FISTA:**

- Uses acceleration for faster convergence
- May temporarily increase objective function
- Best for well-conditioned problems

**Monotone FISTA (MFISTA):**

- Ensures objective function decreases monotonically
- More robust but potentially slower
- Recommended for ill-conditioned problems

**Adaptive Restart (O'Donoghue & Candès, 2015):**

FISTA supports optional adaptive restart to reset the momentum when progress stalls. Set ``adaptive_restart`` to enable:

- ``"function"``: Restart when :math:`f(x^k) > f(x^{k-1})` (objective increased). FISTA only; redundant for MFISTA since it is already monotone.
- ``"gradient"``: Restart when :math:`(y^k - x^{k+1})^T (x^{k+1} - x^k) > 0` (gradient condition). No extra cost; works for both FISTA and MFISTA.

Reference: O'Donoghue, B., & Candès, E. (2015). Adaptive restart for accelerated gradient schemes. Foundations of Computational Mathematics, 15(3), 715-732.

**Example Usage:**

.. code-block:: python

   from pyralysis.optimization.optimizer import FISTA
   from pyralysis.optimization.terms import ChiSquared, L1Norm
   from pyralysis.optimization import ObjectiveFunction

   # Define objective function with L1 sparsity
   term_list = [
       ChiSquared(model_visibility=mv, normalize=True),
       L1Norm(penalization_factor=0.01)  # Single non-differentiable term
   ]
   of = ObjectiveFunction(term_list=term_list, image=image)

   # Create FISTA optimizer
   optim = FISTA(
       objective_function=of,
       parameter=image,
       monotonic=False,  # Use standard FISTA
       adaptive_restart="gradient",  # Optional: O'Donoghue & Candès (2015)
       max_iter=100
   )

   # Run optimization
   reconstructed_image = optim.optimize(verbose=True)

**Key Parameters:**

- ``monotonic``: Whether to use monotone FISTA (default: False)
- ``adaptive_restart``: ``"function"`` or ``"gradient"`` per O'Donoghue & Candès (2015); ``None`` to disable (default)
- ``linesearch``: Line search strategy (default: ``FISTABacktracking`` with
  ``BarzilaiBorweinAdaptiveMin1``; see :doc:`linesearchers` and :doc:`step_size_seeders`)
- ``max_iter``: Maximum iterations
- ``ftol``: Function tolerance for convergence

**Dask note:** FISTA finalizes both the primary iterate and accelerated point each
iteration via :meth:`~pyralysis.optimization.optimizer.Optimizer._finalize_iteration_state`
so graph lineage does not grow with ``max_iter``. See :doc:`optimizers` and
:doc:`performance`.

**Convergence Properties:**

- O(1/k²) convergence rate
- Suitable for large-scale problems
- Requires exactly one non-differentiable term

---

SDMM (Simultaneous Direction Method of Multipliers)
----------------------------------------------------

While FISTA is excellent for single non-differentiable terms, many real-world problems require multiple sparsity constraints or complex regularization schemes. SDMM is designed for problems with multiple non-differentiable terms, making it more flexible than FISTA for complex regularization schemes :cite:`boyd_2011_admm`.

**Mathematical Foundation:**

SDMM solves problems of the form:

.. math::

   \min_x f(x) + \sum_{i=1}^n \lambda_i g_i(x)

where:

- :math:`f(x)` is smooth
- :math:`g_i(x)` are multiple non-differentiable terms
- :math:`\lambda_i` are regularization parameters

**Important:** SDMM requires a differentiable optimizer (such as LBFGS) to handle the smooth term :math:`f(x)`. The differentiable optimizer is passed to SDMM via the ``optimizer`` parameter.

**When to Use SDMM:**

- Multiple sparsity constraints
- Combined regularization (L1 + TV + entropy)
- Complex penalty functions
- Problems where FISTA's single-term limitation is restrictive

**Example Usage:**

.. code-block:: python

   from pyralysis.optimization.optimizer import SDMM, LBFGS
   from pyralysis.optimization.terms import ChiSquared, L1Norm, TSV, Entropy
   from pyralysis.optimization.linesearch import BacktrackingArmijo

   # Define objective function with multiple non-differentiable terms
   term_list = [
       ChiSquared(model_visibility=mv, normalize=True),
       L1Norm(penalization_factor=0.01),      # Sparsity
       TSV(penalization_factor=0.005),        # Smoothness
       Entropy(penalization_factor=0.001)      # Maximum entropy
   ]
   of = ObjectiveFunction(term_list=term_list, image=image)

   # Create differentiable optimizer for smooth terms
   ls = BacktrackingArmijo()
   optim_differentiable = LBFGS(
       parameter=image,
       objective_function=of,
       linesearch=ls,
       max_iter=10
   )

   # Create SDMM optimizer with differentiable optimizer
   optim_sdmm = SDMM(
       parameter=image,
       optimizer=optim_differentiable,  # Required: differentiable optimizer
       max_iter=100,
       objective_function=of,
       rho=0.001
   )

   # Run optimization
   reconstructed_image = optim_sdmm.optimize(verbose=True)

**Key Advantages:**

- Handles multiple non-differentiable terms
- Flexible regularization combinations
- Robust convergence properties
- Suitable for complex imaging problems
- Requires differentiable optimizer for smooth terms (e.g., LBFGS, HagerZhang, etc.)

**Dask note:** SDMM auxiliary variables ``z`` and ``u`` are raw Dask arrays (not
``Parameter`` objects). They are finalized each outer iteration via
:meth:`~pyralysis.optimization.optimizer.Optimizer._finalize_dask_arrays` so graph
lineage does not grow with ``max_iter``. The inner differentiable optimizer still
finalizes the primary iterate via :meth:`~pyralysis.optimization.optimizer.Optimizer._finalize_iteration_state`.
See :doc:`optimizers` and :doc:`performance`.

Regularizers with Proximal Operators
-------------------------------------

Proximal methods (FISTA, SDMM) work with regularizers that have efficient proximal operators. These include both non-differentiable regularizers (which require proximal methods) and smooth regularizers (which can also be used with proximal methods).

**Regularizers with Proximal Operators:**

All regularizers in Pyralysis implement proximal operators, making them compatible with FISTA and SDMM:

- **Non-differentiable** (require proximal methods): ``L1Norm``, ``TVIsotropic``, ``TVAnisotropic``, ``L2Norm``
- **Smooth** (can use gradient-based or proximal methods): ``Huber``, ``Entropy``, ``Tikhonov``, ``TotalFlux``, ``TSV``

**Proximal Operator Characteristics:**

- **L1Norm**: Soft-thresholding (fast, closed-form)
- **TVIsotropic**: Chambolle's algorithm (iterative, rotationally invariant)
- **TVAnisotropic**: Iterative soft-thresholding (faster than isotropic, separable)
- **Huber**: Piecewise linear/quadratic (robust to outliers)
- **L2Norm**: Group soft-thresholding (group sparsity)
- **Entropy**: Exponential scaling (maximum entropy)
- **Tikhonov**: Simple scaling (smoothness)
- **TSV**: FFT-based (edge-preserving smoothness)

.. note::
   For comprehensive documentation on all regularizers, including detailed mathematical formulations, proximal operator implementations, usage examples, and when to use each regularizer, see :doc:`regularization`.

**Example: Combined Regularization**

.. code-block:: python

   from pyralysis.optimization.optimizer import SDMM, LBFGS
   from pyralysis.optimization.terms import ChiSquared, L1Norm, TSV, Huber, TVIsotropic
   from pyralysis.optimization.linesearch import BacktrackingArmijo

   # Multi-term objective function
   term_list = [
       ChiSquared(model_visibility=mv, normalize=True),
       L1Norm(penalization_factor=0.01),           # Sparsity
       TSV(penalization_factor=0.005),              # Smoothness
       # Alternative: Huber(delta=1.0, penalization_factor=0.01),  # Robust to outliers
       # Alternative: TVIsotropic(penalization_factor=0.01)        # True TV (non-smooth)
   ]
   of = ObjectiveFunction(term_list=term_list, image=image)

   # Create differentiable optimizer for smooth terms
   ls = BacktrackingArmijo()
   optim_differentiable = LBFGS(
       parameter=image,
       objective_function=of,
       linesearch=ls,
       max_iter=10
   )

   # Use SDMM for multiple non-differentiable terms
   optim_sdmm = SDMM(
       parameter=image,
       optimizer=optim_differentiable,  # Required: differentiable optimizer
       objective_function=of,
       rho=0.001
   )
   result = optim_sdmm.optimize()

**Algorithm Selection Guide:**

**Use FISTA when:**

- Single non-differentiable term (typically L1 sparsity)
- Large-scale problems requiring fast convergence
- Well-conditioned optimization landscape
- Memory efficiency is important

**Use SDMM when:**

- Multiple non-differentiable terms
- Complex regularization schemes
- Need flexibility in penalty combinations
- FISTA's single-term limitation is restrictive

**Practical Recommendations:**

1. **Start with FISTA** for simple sparsity problems
2. **Upgrade to SDMM** when you need multiple regularizers
3. **Use monotone FISTA** for ill-conditioned problems
4. **Combine L1 + smoothness** for natural image reconstruction

**Performance Considerations:**

**Memory Usage:**

- FISTA: Low memory footprint
- SDMM: Higher memory due to multiple proximal operators

**Convergence Speed:**

- FISTA: O(1/k²) convergence rate
- SDMM: Typically slower but more flexible

**Scalability:**

- Both methods scale well to large images
- FISTA generally faster for single-term problems
- SDMM better for complex regularization

**Proximal Operator Persistence:**

For consistent iteration times in FISTA and SDMM, you can enable proximal operator persistence:

.. code-block:: python

   from pyralysis.optimization.terms.regularizers import L1Norm, Tikhonov

   # Enable persistence for proximal operators
   l1_term = L1Norm(penalization_factor=0.01, persist_proximal=True)
   tikh_term = Tikhonov(penalization_factor=0.001, persist_proximal=True)

   # This maintains consistent iteration times by caching proximal results
   # Similar to persist_gradient for differentiable terms

**Practical Tips:**

- Use appropriate line search parameters
- Monitor convergence with verbose output
- Consider problem conditioning when choosing variants
- Start with moderate regularization parameters

Complete Example: Sparse Reconstruction
---------------------------------------

Complete example of sparse image reconstruction using FISTA:

.. code-block:: python

   import numpy as np
   import dask.array as da
   from pyralysis.optimization.optimizer import FISTA
   from pyralysis.optimization.terms import ChiSquared, L1Norm
   from pyralysis.optimization import ObjectiveFunction
   from pyralysis.reconstruction.parameter import Parameter

   # Create synthetic sparse image
   image_shape = (128, 128)
   true_image = da.zeros(image_shape)
   true_image[32:40, 32:40] = 1.0  # Sparse source
   true_image[80:85, 80:85] = 0.5  # Another sparse source

   # Create parameter
   image = Parameter(data=true_image)

   # Define objective function
   term_list = [
       ChiSquared(model_visibility=mv, normalize=True),
       L1Norm(penalization_factor=0.01)
   ]
   of = ObjectiveFunction(term_list=term_list, image=image)

   # Run FISTA optimization
   optim = FISTA(
       objective_function=of,
       parameter=image,
       monotonic=False,
       max_iter=100,
       ftol=1e-6
   )

   reconstructed_image = optim.optimize(verbose=True)

   print(f"Reconstruction completed")
   print(f"Final image shape: {reconstructed_image.data.shape}")

Conclusion
----------

Compressed sensing provides a powerful framework for reconstructing sparse signals from far fewer measurements than traditional methods require. By combining the theoretical foundations (sparsity, incoherence, RIP) with practical optimization tools (proximal operators, FISTA, SDMM), we can solve challenging radio interferometric imaging problems efficiently.

The key insights are:

- **Sparsity** enables reconstruction from fewer measurements
- **L1 relaxation** makes the problem tractable
- **Proximal operators** handle non-differentiability
- **FISTA/SDMM** provide efficient optimization algorithms

For radio interferometry, compressed sensing is particularly valuable when complete :math:`(u,v)` coverage is impossible, making it an essential tool for modern radio astronomy.

References
----------

.. bibliography::
   :filter: docname in docnames

See Also
--------
- :doc:`regularization` - Comprehensive guide to all regularizers and their properties
- :doc:`optimizers` - Gradient-based methods (for smooth regularizers)
- :doc:`linesearchers` - ``FISTABacktracking`` and Armijo line search
- :doc:`step_size_seeders` - Proximal and gradient BB seeders
- :doc:`optimization` - Workflow overview
- :doc:`inverse_problem` - Mathematical foundations of inverse problems
- :doc:`autoapi/index` - API reference

----

:doc:`regularization` | :doc:`optimization` | :doc:`inverse_problem`