Compressed Sensing & Proximal Methods

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 [Candès et al., 2006, Donoho, 2006]. In radio interferometric imaging, this is particularly powerful because astronomical sources are often sparse in appropriate bases [Carrillo et al., 2014, Wiaux et al., 2009].

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 \(B\), this requires at least \(2B\) samples per second. In radio interferometry, this translates to needing complete coverage of the \((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 [Candès et al., 2006, Donoho, 2006]. Instead of needing \(n\) measurements for an \(n\)-dimensional signal, we only need roughly \(k \log(n/k)\) measurements, where \(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 \((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:

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

where:

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

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

Formulating the Compressed Sensing Problem:

In traditional linear inverse problems, we solve:

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

where \(\mathbf{A}\) is the measurement matrix and \(\mathbf{b}\) are the measurements.

For compressed sensing, this becomes:

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

where:

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

This formulation allows us to work directly with the sparse coefficients \(\mathbf{x}\) rather than the full signal \(\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):

\[\|\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:

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

This relaxation has several key advantages [Donoho, 2006]:

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:

\[\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
\(\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 \(\mathbf{A}\) must be incoherent with respect to the sparsity basis \(\mathbf{\Phi}\). Incoherence means that the measurement vectors cannot be too aligned with the sparsity basis vectors.

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

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

where \(\mathbf{a}_i\) are the rows of \(\mathbf{A}\) and \(\mathbf{\phi}_j\) are the columns of \(\mathbf{\Phi}\).

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

2. Sparsity Level Requirement: The number of measurements \(p\) must satisfy:

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

where:

\(k\) is the number of non-zero elements in \(\mathbf{x}\)
\(n\) is the dimension of \(\mathbf{x}\)

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

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

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

\[(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 \(\delta_{2k} < 0.414\).

Application to Radio Interferometry:

In radio interferometry, these conditions translate to:

\(\mathbf{A}\) represents the sampling pattern in the \((u,v)\) plane

\(\mathbf{\Phi}\) could be the Fourier basis (for point sources) or wavelet basis (for extended sources)

\(\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 \((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 \(g\) at point \(\mathbf{x}\) is defined as:

\[\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 \(\lambda > 0\) is a parameter controlling the strength of the penalty.

Key Properties:

Nonexpansive: \(\|\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 \(g(\mathbf{x}) = \|\mathbf{x}\|_1\), the proximal operator is:

\[\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 \(g(\mathbf{x}) = \frac{1}{2}\|\mathbf{x}\|_2^2\), the proximal operator is:

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

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

\[\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 \(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 \(k\) is the number of iterations [Beck and Teboulle, 2009]. It’s particularly effective for problems with exactly one non-differentiable term.

Mathematical Foundation:

FISTA solves problems of the form:

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

where:

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

Algorithm Steps:

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

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

\[\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 \(f(x^k) > f(x^{k-1})\) (objective increased). FISTA only; redundant for MFISTA since it is already monotone.
"gradient": Restart when \((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:

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 Line Searchers and 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 _finalize_iteration_state() so graph lineage does not grow with max_iter. See Optimizers and Performance & Scalability.

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 [Boyd et al., 2011].

Mathematical Foundation:

SDMM solves problems of the form:

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

where:

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

Important: SDMM requires a differentiable optimizer (such as LBFGS) to handle the smooth term \(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:

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 _finalize_dask_arrays() so graph lineage does not grow with max_iter. The inner differentiable optimizer still finalizes the primary iterate via _finalize_iteration_state(). See Optimizers and Performance & Scalability.

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 Regularization Methods.

Example: Combined Regularization

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:

Start with FISTA for simple sparsity problems
Upgrade to SDMM when you need multiple regularizers
Use monotone FISTA for ill-conditioned problems
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:

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:

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 \((u,v)\) coverage is impossible, making it an essential tool for modern radio astronomy.

References#

[1]

Amir Beck and Marc Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183–202, 2009. doi:10.1137/080716542.

[2]

Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning, 3(1):1–122, 2011. doi:10.1561/2200000016.

[3] (1,2)

Emmanuel J. Candès, Justin Romberg, and Terence Tao. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory, 52(2):489–509, 2006. doi:10.1109/TIT.2005.862083.

[4]

Rafael E. Carrillo, Jason D. McEwen, and Yves Wiaux. Purify: a new approach to radio-interferometric imaging. Monthly Notices of the Royal Astronomical Society, 439(4):3591–3604, 2014. doi:10.1093/mnras/stu202.

[5] (1,2,3)

David L. Donoho. Compressed sensing. IEEE Transactions on Information Theory, 52(4):1289–1306, 2006. doi:10.1109/TIT.2006.871582.

[6]

Yves Wiaux, Laurent Jacques, Gilles Puy, Anna M. M. Scaife, and Pierre Vandergheynst. Compressed sensing imaging techniques for radio interferometry. Monthly Notices of the Royal Astronomical Society, 395(3):1733–1742, 2009. doi:10.1111/j.1365-2966.2009.14665.x.

Compressed Sensing & Proximal Methods

Contents

Compressed Sensing & Proximal Methods#

Introduction#

Mathematical Foundation of Compressed Sensing#

From Theory to Optimization#

Theoretical Guarantees#

Proximal Operators: The Key to Practical Implementation#

Proximal Optimization Algorithms#

FISTA (Fast Iterative Shrinkage-Thresholding Algorithm)#

SDMM (Simultaneous Direction Method of Multipliers)#

Regularizers with Proximal Operators#

Complete Example: Sparse Reconstruction#

Conclusion#

References#

See Also#