Plug-and-Play Volumetric Reconstruction for Compressive Sensing Light-Sheet Microscopy5
July 02, 2026
Abstract
We investigate volumetric reconstruction for compressive sensing light-sheet microscopy (CS-LSM), where fast volumetric imaging is achieved by encoding multiple axial planes into each camera exposure. To recover the underlying volume from highly multiplexed measurements, we propose a plug-and-play (PnP) framework that flexibly incorporates any user-specified denoiser into the reconstruction process. Building on a slice-based formulation, we further introduce an axial-coupled model that exploits correlations between adjacent slices to improve volumetric continuity. For efficient computation, we derive a Woodbury-based update for the data-consistency step in both the slice-based and axial-coupled formulations, and employ a Gauss–Seidel sweep for the denoising step in the axial-coupled model. Under a weakly convex regularization assumption, we establish subsequential convergence of the proposed algorithm. Experiments on synthetic and real zebrafish-heart data demonstrate that the proposed framework successfully recovers cellular structures from compressed measurements, and provide practical insights into the comparative performance of commonly used denoisers within the PnP framework under the CS-LSM setup.
Compressive Sensing, Light-Sheet Microscopy, Plug-and-Play, Image Reconstruction
68U10, 65K10, 65F22, 94A08, 92C55
1 Introduction↩︎
Biological processes such as cardiac contraction and intracardiac flow evolve in three-dimensional (3D) space and over extremely short time scales [1]. Capturing these dynamics in vivo requires volumetric imaging methods that achieve cellular resolution over organ-scale fields of view while maintaining sufficiently high temporal resolution [2]. At the same time, light exposure must be carefully controlled to limit photobleaching and phototoxicity [3], as these measurements are often performed in living specimens.
In fluorescence microscopy, volumetric imaging speed, spatial resolution, field of view, and excitation burden are constrained by fundamental trade-offs, thereby hindering their simultaneous achievement. Existing 3D imaging methods for beating hearts include confocal laser scanning microscopy (CLSM) [4]–[6], two-photon microscopy (TPM) [7]–[9], and light-sheet microscopy (LSM) [10]–[15]. CLSM and TPM offer strong optical sectioning, but their point- or line-scanning acquisition limits the rate at which full volumes can be recorded and may lead to substantial photobleaching, phototoxicity, or thermal load during prolonged high-speed imaging [16], [17]. In contrast, LSM illuminates only the imaging plane and detects fluorescence from an orthogonal direction, positioning it as an attractive approach for fast volumetric imaging with reduced light dose. Despite these advantages, conventional LSM struggles to capture dynamics that evolve on millisecond timescales across 3D volumes, and retrospective synchronization or prospective optical gating is often employed to compensate [12]–[15]. Approaches such as rapid beam steering [18] or multiplane recording [19] partially mitigate this limitation but typically rely on specialized optical instrumentation and detection hardware. More fundamentally, detector bandwidth imposes hard limits on volumetric imaging performance, as higher acquisition rates come at the expense of spatial sampling, field of view, or signal quality.
This trade-off motivates coded and compressive acquisition strategies, and has recently driven the integration of LSM with compressive sensing (CS) [20], [21] for more efficient volumetric encoding. While such coded and compressive volumetric imaging approaches show clear promise for improving volumetric imaging speed and reducing the number of required measurements, achieving fast in vivo volumetric imaging at cellular resolution remains challenging. For instance, a spatially modulated LSM recover volumes from patterned illumination using CS, but it relies on multiple sequential acquisitions, limiting volumetric speed [22]. Snapshot temporal compressive LSM can recover multiple frames from a single measurement, but it is typically designed for a single light-sheet plane rather than a full 3D volume [23]. In addition, light field microscopy enables single-shot volumetric capture by encoding angular information onto a single sensor; however, the associated trade-off between angular and spatial sampling often reduces effective resolution, introduces on-focal artifacts, or requires customized optical components and extensive training data for neural network-based reconstruction, making cellular-scale in vivo imaging challenging [24]–[27]. Taken together, these studies indicate that fast acquisition alone is insufficient; when measurements are strongly multiplexed and undersampled, effective reconstruction becomes a central component of the imaging pipeline.
Addressing the acquisition side of this pipeline, prior work by our team developed a CS-LSM platform for compressed light-sheet acquisition [28]. This paper focuses on the reconstruction problem; for details on hardware design, calibration, and acquisition protocol, we refer the reader to [29]. In the CS-LSM system, axial scanning is synchronized with spatial light modulation via a digital micromirror device (DMD). Following the incoherent measurement framework of CS [21], random binary masks are used to encode fluorescence from multiple depth planes into a single exposure, enabling volumetric imaging at rates of 200 volumes per second while maintaining cellular resolution [28]. It also reduces the number of recorded measurements and the associated storage requirements relative to non-compressed acquisition, thereby improving data efficiency in high-speed volumetric imaging [29].
While this hardware platform enables high-speed compressed acquisition, the resulting measurements are highly multiplexed and undersampled, necessitating robust and efficient reconstruction algorithms to recover high-quality volumetric reconstructions. The main objective of this paper is therefore to develop a flexible reconstruction framework tailored to this hardware system. Although motivated by the present CS-LSM platform, the resulting methodology may be applicable to a broader class of coded, multiplexed, and computational volumetric imaging settings.
To this end, we build on a plug-and-play (PnP) framework [30] and solve the model using the alternating direction method of multipliers (ADMM) [31], which we refer to as PnP-ADMM. We first formulate a slice-based reconstruction model, in which each slice is reconstructed independently. This design exploits the diagonal structure of the DMD masks, enabling efficient data-consistency updates via Woodbury-based inversion. To better preserve volumetric continuity, we subsequently introduce an axially coupled reconstruction model that enforces smoothness between adjacent slices along the \(z\) direction. This axial coupling is designed to leverage inter-slice correlations within the reconstructed volume, rather than temporal correlations across different cardiac phases. In contrast to full 3D regularization or volumetric patch-based reconstruction methods [32]–[34], we incorporate axial coupling within the denoising step using a Gauss–Seidel sweep. This approach is computationally efficient and, importantly, enables us to establish convergence guarantees for the proposed algorithm.
The name plug-and-play reflects the fact that one subproblem is formulated as a denoising step, allowing a wide range of denoisers to be incorporated without altering the overall reconstruction framework. This flexibility has also been leveraged in hyperspectral unmixing [35], limited-angle tomography [36], and autonomous driving [37]. Here, we investigate the performance of classical denoisers, such as Tikhonov [38], total variation (TV) [39]–[42], and block-matching and 3D filtering (BM3D) [43], as well as deep learning (DL)-based denoisers, including the denoising convolutional neural network (DnCNN) [44], the fast and flexible denoising network (FFDNet) [45], and the deep residual U-Net denoiser (DRUNet) [46]. Experiments on both synthetic and real zebrafish-heart data demonstrate the effectiveness of the proposed framework and provide practical insights into the strengths and limitations of these off‑the‑shelf denoisers.
For convex regularizers such as Tikhonov and TV, the convergence of PnP-ADMM follows directly from classical ADMM theories [31], [47]. However, this guarantee does not readily extend to more advanced denoisers, such as BM3D and DL-based methods, for which an explicit underlying prior is typically unavailable, let alone one that can be verified to be convex. To address this gap, Chan et al. [48] established convergence of PnP algorithms under the assumption that the denoiser is nonexpansive, while Ryu et al. [49] relaxed this requirement by imposing nonexpansiveness on the residual operator, at the cost of requiring strong convexity of the data-fidelity term, which is not satisfied in the CS-LSM setting. In this work, we establish subsequential convergence of the PnP-ADMM algorithm for the proposed axial-coupled model for a weakly convex regularization term.
The main contributions of this paper are summarized as follows:
We introduce an axial-coupled reconstruction model that leverages inter-slice correlations along the \(z\) direction, leading to improved reconstruction performance compared with slice-based model.
We tailor the PnP-ADMM framework to the CS-LSM setup through efficient algorithmic updates, including a Woodbury-based reformulation of the linear subproblem and a Gauss–Seidel update strategy for the axial-coupled model.
We analyze the convergence properties of the proposed algorithm, establishing subsequential convergence under the assumption that the image prior is weakly convex.
Experiments on synthetic and real data provide practical guidance for selecting off-the-shelf denoisers in the CS-LSM reconstruction problem.
The rest of this paper is organized as follows. Section 2 introduces the CS-LSM forward model and the proposed PnP-ADMM reconstruction framework, including both slice-based and axial-coupled models. Section 3 discusses convergence properties of the proposed framework. Section 4 presents numerical results on synthetic and real datasets. Finally, Section 5 concludes the paper and outlines future work.
2 Proposed Approaches↩︎
In this section, we first introduce the CS-LSM forward model in Section 2.1. We then present a slice-based reconstruction model in Section 2.2, where a 3D volume is recovered through a collection of 2D subproblems. Next, in Section 2.3, we extend this formulation to an axial-coupled model by enforcing inter-slice smoothness along the \(z\) direction. Both models are solved within an ADMM framework. Notably, the regularization subproblem is interpreted as a denoising step, reflecting the plug-and-play nature of the framework.
2.1 Forward model and notation↩︎
We adopt the CS-LSM acquisition model from [28] and focus on the reconstruction problem. Specifically, we obtain compressed image data through a set of binary masks, as illustrated in Fig. 1. The top row of Fig. 1 shows a sequence of 3D cardiac volumes acquired at \(M\) distinct time points during one heartbeat, indexed by \(t_1,t_2,\dots,t_M\). In this paper, we focus on reconstruction at a fixed time point, recovering one single 3D volume at a time. Extending this across time yields a 4D sequence, while explicit temporal coupling is left for future work. At a time point, the underlying 3D volume is acquired through \(N\) compressed camera shots over one scanning period. As shown in Fig. 1, each shot corresponds to one axial block of the volume, and each axial block contains \(R\) axial slices along the \(z\) direction, where \(R\) is the compression ratio, namely, the number of axial slices encoded within a single camera exposure. Consequently, the full reconstructed volume contains \(S:=NR\) axial slices in total. We denote these unknown slices by \(\{\boldsymbol{v}_n\}_{n=1}^{S}\), where each \(\boldsymbol{v}_n\in\mathbb{R}^{p}\) represents one vectorized 2D slice6 and \(p\) denotes the total number of pixels in each 2D slice.
In the CS acquisition process, the DMD applies \(R\) binary masks sequentially during each camera exposure. Within the \(j\)th shot (\(j=1, \cdots, N\)), the \(R\) axial slices in the corresponding section group are modulated by distinct binary masks through element-wise (Hadamard) products and then summed on the detector, producing a single compressed 2D measurement, which we vectorize as \(\boldsymbol{b}_j\in\mathbb{R}^{p}\). Let \(\boldsymbol{\phi}_r\in\mathbb{R}^{p\times p}\) denote the diagonal masking operator for the \(r\)th slice (\(r=1, \cdots, R)\), so that the Hadamard modulation can be written as a matrix-vector product. Then the \(j\)th compressed measurement can be mathematically modeled as \[\label{eq:forward95model} \boldsymbol{b}_j=\sum_{r=1}^{R}\boldsymbol{\phi}_r \boldsymbol{v}_{(j-1)R+r}+\boldsymbol{\eta}_j,\qquad j=1,\dots,N,\tag{1}\] where \(\boldsymbol{\eta}_j\) denotes measurement noise. Thus, the reconstruction task is to recover the full stack of \(S\) axial slices \(\{\boldsymbol{v}_n\}_{n=1}^{S}\) from only \(N\) compressed measurements stored in \(\{\boldsymbol{b}_j\}_{j=1}^{N}\), which is a highly underdetermined inverse problem. Moreover, since adjacent axial slices in the volume are often strongly correlated, it is beneficial to exploit such inter-slice structure during reconstruction, which motivates the axial-coupled model introduced in Section 2.3.
2.2 Slice-based reconstruction↩︎
We begin with a slice-based model, in which each axial slice is regularized individually. Assuming that the measurement noise can be well approximated as additive Gaussian7, we employ a least-squares data-fidelity term and consider the following optimization problem, \[\label{eq:model-slice} \min_{\{\boldsymbol{v}_n\}_{n=1}^{S}} \frac{1}{2}\sum_{j=1}^{N}\left\|\boldsymbol{b}_j-\sum_{r=1}^{R}\boldsymbol{\phi}_r \boldsymbol{v}_{(j-1)R+r}\right\|^2 +\lambda\sum_{n=1}^{S}\psi(\boldsymbol{v}_n),\tag{2}\] where \(\psi(\cdot)\) denotes a slice-wise regularization term and \(\lambda>0\) is a weighting parameter.
To solve 2 , we introduce auxiliary variables \(\{\boldsymbol{u}_n\}_{n=1}^{S}\) and rewrite the model as \[\label{eq:model-slice-split} \begin{align} \min_{\{\boldsymbol{v}_n,\boldsymbol{u}_n\}} \quad & \frac{1}{2}\sum_{j=1}^{N}\left\|\boldsymbol{b}_j-\sum_{r=1}^{R}\boldsymbol{\phi}_r \boldsymbol{v}_{(j-1)R+r}\right\|^2 +\lambda\sum_{n=1}^{S}\psi(\boldsymbol{u}_n) \\ \text{subject to}\quad & \boldsymbol{v}_n=\boldsymbol{u}_n,\qquad n=1,\dots,S. \end{align}\tag{3}\] The corresponding augmented Lagrangian is \[\begin{align} \mathcal{L}_{\rho}(\{\boldsymbol{u}_n\},\{\boldsymbol{v}_n\},\{\boldsymbol{d}_n\}) :=&\; \frac{1}{2}\sum_{j=1}^{N}\left\|\boldsymbol{b}_j-\sum_{r=1}^{R}\boldsymbol{\phi}_r \boldsymbol{v}_{(j-1)R+r}\right\|^2\notag \\ &+\sum_{n=1}^{S}\left( \lambda\psi(\boldsymbol{u}_n) +\frac{\rho}{2}\|\boldsymbol{v}_n-\boldsymbol{u}_n+\boldsymbol{d}_n\|^2 -\frac{\rho}{2}\|\boldsymbol{d}_n\|^2 \right), \end{align}\] where \(\{\boldsymbol{d}_n\}_{n=1}^{S}\) are scaled dual variables and \(\rho>0\) is the penalty parameter. The ADMM iterations are given by \[\label{eq:ADMM} \left\{ \begin{align} \{\boldsymbol{u}_n^{k+1}\} &= \arg\min_{\{\boldsymbol{u}_n\}} \mathcal{L}_{\rho}(\{\boldsymbol{u}_n\},\{\boldsymbol{v}_n^k\},\{\boldsymbol{d}_n^k\}),\\ \{\boldsymbol{v}_n^{k+1}\} &= \arg\min_{\{\boldsymbol{v}_n\}} \mathcal{L}_{\rho}(\{\boldsymbol{u}_n^{k+1}\},\{\boldsymbol{v}_n\},\{\boldsymbol{d}_n^k\}),\\ \boldsymbol{d}_n^{k+1} &= \boldsymbol{d}_n^k+\boldsymbol{v}_n^{k+1}-\boldsymbol{u}_n^{k+1}, \end{align} \right.\tag{4}\] where \(n=1,\dots,S\) and the superscript \(k\) counts the iteration number.
The \(\boldsymbol{u}\)-subproblem decouples across slices: \[\label{eq:u-sub} \boldsymbol{u}_n^{k+1} = \arg\min_{\boldsymbol{z}}\; \lambda\psi(\boldsymbol{z})+\frac{\rho}{2}\|\boldsymbol{z}-(\boldsymbol{v}_n^{k}+\boldsymbol{d}_n^k)\|^2, \qquad n=1,\dots,S.\tag{5}\] From the PnP perspective, this subproblem is interpreted as a denoising step applied to the input \(\boldsymbol{v}_n^{k}+\boldsymbol{d}_n^k\). Many existing denoising approaches, or denoisers for short, can be formulated as \[\label{eq:denoiser} \mathcal{K}_\sigma(\boldsymbol{f}) := \arg\min_{\boldsymbol{z}} \psi(\boldsymbol{z})+ \frac{1}{2\sigma^2}\|\boldsymbol{z}-\boldsymbol{f}\|^2,\tag{6}\] where \(\sigma>0\) is an input noise-level parameter common to methods such as BM3D, DnCNN, FFDNet, and DRUNet, and typically corresponds to the standard deviation of the noise. Dividing 5 by \(\rho\) and comparing with the denoising formulation 6 , we identify the denoiser noise level as \(\sigma_\lambda=\sqrt{\lambda/\rho}\) and adopt the following update: \[\label{eq:u-pnp-slice} \boldsymbol{u}_n^{k+1}=\mathcal{K}_{\sigma_\lambda}(\boldsymbol{v}_n^{k}+\boldsymbol{d}_n^k),\tag{7}\] where \(\mathcal{K}_{\sigma_\lambda}\) denotes the selected denoiser with its internal noise level \(\sigma_\lambda\) controlled by \(\lambda.\) When \(\mathcal{K}_{\sigma_\lambda}\) is induced by an explicit regularizer, 5 reduces to a proximal update. As an example, Tikhonov regularization with \(\psi(\boldsymbol{z})=\|\boldsymbol{z}\|^2\) yields \[\label{eq:tikhonov-slice} \boldsymbol{u}_n^{k+1} = \frac{\rho(\boldsymbol{v}_n^{k}+\boldsymbol{d}_n^k)}{2\lambda+\rho}.\tag{8}\]
For the \(\boldsymbol{v}\)-subproblem, the variables associated with each compressed shot can be updated shotwise. Fix \(j\in\{1,\dots,N\}\) and define \[\boldsymbol{c}_r^{(j),k+1} := \boldsymbol{u}_{(j-1)R+r}^{k+1} - \boldsymbol{d}_{(j-1)R+r}^{k}, \qquad r=1,\dots,R.\] Let \[\boldsymbol{v}^{(j)} := [\boldsymbol{v}_{(j-1)R+1}^\top,\dots,\boldsymbol{v}_{jR}^\top]^\top \in\mathbb{R}^{pR}, \qquad \boldsymbol{c}^{(j),k+1} := [(\boldsymbol{c}_1^{(j),k+1})^\top,\dots,(\boldsymbol{c}_R^{(j),k+1})^\top]^\top \in\mathbb{R}^{pR},\] and \[\Phi:=[\boldsymbol{\phi}_1,\dots,\boldsymbol{\phi}_R]\in\mathbb{R}^{p\times pR}.\] Then the corresponding subproblem becomes \[\label{eq:v-sub-1} \min_{\boldsymbol{v}^{(j)}}\; \frac{1}{2}\|\boldsymbol{b}_j-\Phi \boldsymbol{v}^{(j)}\|^2+\frac{\rho}{2}\|\boldsymbol{v}^{(j)}-\boldsymbol{c}^{(j),k+1}\|^2.\tag{9}\] Its closed-form solution is \[\label{eq:v-update} \boldsymbol{v}^{(j),k+1} = (\Phi^{\top}\Phi+\rho I_{pR})^{-1} (\Phi^{\top}\boldsymbol{b}_j+\rho\boldsymbol{c}^{(j),k+1}).\tag{10}\]
The matrix inverse in 10 can be simplified by the Woodbury matrix identity \[\label{Woodbury} (I+UV)^{-1}=I-U(I+VU)^{-1}V.\tag{11}\] Applying 11 yields the equivalent update \[\label{eq:v-update-fast} \boldsymbol{v}^{(j),k+1} = \boldsymbol{c}^{(j),k+1} +\Phi^{\top}(\rho I_p+\Phi\Phi^{\top})^{-1} (\boldsymbol{b}_j-\Phi\boldsymbol{c}^{(j),k+1}).\tag{12}\] Equation 12 is computationally attractive because it replaces the inversion of the \(pR\times pR\) matrix in 10 by the inversion of the \(p\times p\) matrix \(\rho I_p+\Phi\Phi^{\top}\), which has a simple structure. Indeed, since each \(\boldsymbol{\phi}_r\) is diagonal, then the product \(\boldsymbol{\phi}_r\boldsymbol{\phi}_r^{\top}\) is diagonal and \(\Phi\Phi^{\top}=\sum_{r=1}^R \boldsymbol{\phi}_r\boldsymbol{\phi}_r^{\top}\) is diagonal as well. Therefore, \((\rho I_p+\Phi\Phi^{\top})^{-1}\) can be computed element-wise. More specifically, if \(\phi_r(i)\) denotes the \(i\)th diagonal entry of \(\boldsymbol{\phi}_r\), then \[\label{eq:phi-diag} [\Phi\Phi^{\top}]_{ii}=\sum_{r=1}^R \phi_r(i)^2.\tag{13}\] For binary masks, \(\phi_r(i)^2=\phi_r(i)\), so the quantity in 13 is simply the number of masks that are active at pixel \(i\). Therefore, the data-consistency step in 12 does not require any matrix inversion and can instead be implemented using inexpensive pointwise operations. Specifically, for a fixed shot \(j\), 12 admits a pixel-wise interpretation: for each detector pixel location \(i\), define the \(R\)-dimensional vectors \[\boldsymbol{m}(i) := [\phi_1(i),\dots,\phi_R(i)]^{\top}\in\mathbb{R}^R, \qquad \widehat{\boldsymbol{c}}^{(j),k+1}(i) := [c_1^{(j),k+1}(i),\dots,c_R^{(j),k+1}(i)]^{\top}\in\mathbb{R}^R,\] and let \[\widehat{\boldsymbol{v}}^{(j),k+1}(i) := [v_{(j-1)R+1}^{k+1}(i),\dots,v_{jR}^{k+1}(i)]^{\top}\in\mathbb{R}^R,\] where \(\widehat{\boldsymbol{c}}^{(j),k+1}(i)\) and \(\widehat{\boldsymbol{v}}^{(j),k+1}(i)\) collect the current and updated values of the \(R\) slices at pixel location \(i\) within the \(j\)th shot, respectively. Then 12 can be written as \[\label{eq:v-update-pixel} \widehat{\boldsymbol{v}}^{(j),k+1}(i) = \widehat{\boldsymbol{c}}^{(j),k+1}(i) + \boldsymbol{m}(i) \frac{ b_j(i)-\boldsymbol{m}(i)^{\top}\widehat{\boldsymbol{c}}^{(j),k+1}(i) }{ \rho+\|\boldsymbol{m}(i)\|_2^2 }, \quad i=1, \cdots, p.\tag{14}\] This expression shows that the data-consistency step decouples across spatial pixels: different detector pixels can be updated independently, while the coupling occurs only among the \(R\) slices sharing the same detector location. In other words, 14 is a residual correction along the local mask direction \(\boldsymbol{m}(i)\), scaled by the normalization factor \(\rho+\|\boldsymbol{m}(i)\|_2^2\).
2.3 Axial-coupled reconstruction↩︎
While the slice-based model treats each axial slice independently, adjacent slices in a 3D volume are often strongly correlated. Ignoring this inter-slice structure may lead to discontinuities or inconsistent reconstructions along the \(z\) direction. A natural remedy is to apply full 3D regularization directly to the reconstructed volume, such as 3D TV [32], BM4D [33] (a volumetric extension of BM3D), or temporal non-local means [34]. We incorporate a quadratic axial coupling term that penalizes differences between neighboring slices. This pairwise coupling admits an efficient Gauss–Seidel (GS) update strategy for the denoising step, for which we establish subsequential convergence in Section 3. Concretely, the axial-coupled model reads \[\label{eq:model-axial} \min_{\{\boldsymbol{v}_n\}_{n=1}^{S}} \frac{1}{2}\sum_{j=1}^{N}\left\|\boldsymbol{b}_j-\sum_{r=1}^{R}\boldsymbol{\phi}_r \boldsymbol{v}_{(j-1)R+r}\right\|^2 +\lambda\sum_{n=1}^{S}\psi(\boldsymbol{v}_n) +\frac{\gamma}{2}\sum_{n=1}^{S}\|\boldsymbol{v}_n-\boldsymbol{v}_{n-1}\|^2,\tag{15}\] where \(\gamma>0\) controls the strength of axial coupling. For notational convenience, we impose periodic boundary conditions, \(\boldsymbol{v}_0=\boldsymbol{v}_S\).
As in the slice-based case, we introduce auxiliary variables \(\{\boldsymbol{u}_n\}_{n=1}^{S}\) and rewrite 15 as \[\label{eq:model-axial-split} \begin{align} \min_{\{\boldsymbol{v}_n,\boldsymbol{u}_n\}} \quad & \frac{1}{2}\sum_{j=1}^{N}\left\|\boldsymbol{b}_j-\sum_{r=1}^{R}\boldsymbol{\phi}_r \boldsymbol{v}_{(j-1)R+r}\right\|^2 +\lambda\sum_{n=1}^{S}\psi(\boldsymbol{u}_n) +\frac{\gamma}{2}\sum_{n=1}^{S}\|\boldsymbol{u}_n-\boldsymbol{u}_{n-1}\|^2 \\ \text{subject to}\quad & \boldsymbol{v}_n=\boldsymbol{u}_n,\qquad n=1,\dots,S. \end{align}\tag{16}\] The corresponding augmented Lagrangian is \[\begin{align} \mathcal{L}_{\rho}(\{\boldsymbol{u}_n\},\{\boldsymbol{v}_n\},\{\boldsymbol{d}_n\}) :=&\; \frac{1}{2}\sum_{j=1}^{N}\left\|\boldsymbol{b}_j-\sum_{r=1}^{R}\boldsymbol{\phi}_r \boldsymbol{v}_{(j-1)R+r}\right\|^2 \notag\\ +\sum_{n=1}^{S}&\left( \lambda\psi(\boldsymbol{u}_n) +\frac{\rho}{2}\|\boldsymbol{v}_n-\boldsymbol{u}_n+\boldsymbol{d}_n\|^2 -\frac{\rho}{2}\|\boldsymbol{d}_n\|^2 \right) +\frac{\gamma}{2}\sum_{n=1}^{S}\|\boldsymbol{u}_n-\boldsymbol{u}_{n-1}\|^2, \label{eq:aug-axial} \end{align}\tag{17}\] where \(\{\boldsymbol{d}_n\}_{n=1}^{S}\) are scaled dual variables and \(\rho>0\) is the penalty parameter.
Following the ADMM scheme 4 , we first update the auxiliary variables \(\{\boldsymbol{u}_n\}_{n=1}^{S}\) with \(\{\boldsymbol{v}_n^k\}\) and \(\{\boldsymbol{d}_n^k\}\) fixed. Due to the axial coupling term, the \(\boldsymbol{u}\)-subproblem no longer decouples across slices: \[\label{eq:u-sub-global-axial} \{\boldsymbol{u}^{k+1}_n\}_{n=1}^{S}=\arg\min_{\{\boldsymbol{z}_n\}_{n=1}^{S}} \lambda\sum_{n=1}^{S}\psi(\boldsymbol{z}_n) +\frac{\rho}{2}\sum_{n=1}^{S}\|\boldsymbol{z}_n-(\boldsymbol{v}_n^{k}+\boldsymbol{d}_n^k)\|^2 +\frac{\gamma}{2}\sum_{n=1}^{S}\|\boldsymbol{z}_n-\boldsymbol{z}_{n-1}\|^2,\tag{18}\] that is, each \(\boldsymbol{u}_n\) depends on its neighboring slices \(\boldsymbol{u}_{n-1}\) and \(\boldsymbol{u}_{n+1}.\) We adopt a GS strategy, updating one slice at a time in a cyclic sweep with periodic boundary conditions, i.e., \(\boldsymbol{u}_0=\boldsymbol{u}_S\) and \(\boldsymbol{u}_{S+1}=\boldsymbol{u}_1\), where each slice is updated with its two neighbors fixed at their most recently computed values. More precisely, for \(n=1,\dots,S\), define \[\label{eq:u-periodic} \boldsymbol{u}_{n,-}^{k} := \begin{cases} \boldsymbol{u}_{n-1}^{k+1}, & n=2,\dots,S,\\ \boldsymbol{u}_S^{k}, & n=1, \end{cases} \qquad \boldsymbol{u}_{n,+}^{k} := \begin{cases} \boldsymbol{u}_{n+1}^{k}, & n=1,\dots,S-1,\\ \boldsymbol{u}_1^{k+1}, & n=S. \end{cases}\tag{19}\] Then each slice is updated by solving \[\boldsymbol{u}_n^{k+1} = \arg\min_{\boldsymbol{z}}\; \lambda\psi(\boldsymbol{z}) +\frac{\rho}{2}\|\boldsymbol{z}-(\boldsymbol{v}_n^{k}+\boldsymbol{d}_n^k)\|^2 +\frac{\gamma}{2}\|\boldsymbol{z}-\boldsymbol{u}_{n,-}^{k}\|^2 +\frac{\gamma}{2}\|\boldsymbol{z}-\boldsymbol{u}_{n,+}^{k}\|^2 . \label{eq:axial95u95sub}\tag{20}\]
Completing the square and defining \[\label{eq:g-update} \boldsymbol{g}_n^k = \frac{\rho(\boldsymbol{v}_n^{k}+\boldsymbol{d}_n^k)+\gamma(\boldsymbol{u}_{n,-}^{k}+\boldsymbol{u}_{n,+}^{k})}{\rho+2\gamma}\in\mathbb{R}^p,\tag{21}\] the subproblem 20 reduces to \[\label{eq:u-update-axial} \boldsymbol{u}_n^{k+1} = \arg\min_{\boldsymbol{z}}\; \lambda\psi(\boldsymbol{z})+\frac{\rho+2\gamma}{2}\|\boldsymbol{z}-\boldsymbol{g}_n^k\|^2.\tag{22}\] From the PnP perspective, 22 can be interpreted as a denoising step applied to the axial-coupled quantity \(\boldsymbol{g}_n^k\). Normalizing 22 by \(\rho+2\gamma\) and comparing with 6 leads to the denoising update: \[\label{eq:u-pnp-axial} \boldsymbol{u}_n^{k+1} = \mathcal{K}_{\sigma_\lambda}(\boldsymbol{g}_n^k),\tag{23}\] with \(\sigma_\lambda=\sqrt{\frac{\lambda}{\rho+2\gamma}}\). When \(\mathcal{K}_{\sigma_\lambda}\) is induced by an explicit regularizer, for example, \(\psi(\boldsymbol{z})=\|\boldsymbol{z}\|^2\), the Tikhonov update can be written as the closed-form update \[\label{eq:tikhonov-axial} \boldsymbol{u}_n^{k+1} = \frac{ \rho(\boldsymbol{v}_n^{k}+\boldsymbol{d}_n^k) +\gamma(\boldsymbol{u}_{n,-}^{k}+\boldsymbol{u}_{n,+}^{k}) }{2\lambda+\rho+2\gamma}.\tag{24}\]
With \(\{\boldsymbol{u}_n^{k+1}\}_{n=1}^{S}\) fixed, the \(\boldsymbol{v}\)- and \(\boldsymbol{d}\)-updates are the same as in the slice-based model. We summarize the axial-coupled PnP-ADMM scheme in Algorithm 2. When \(\gamma=0\), the algorithm reduces to the slice-based model.
3 Convergence Analysis↩︎
In this section, we study convergence properties of the proposed Algorithm 2. To facilitate the analysis, we recast the summation-based model 15 into an equivalent matrix-vector formulation and introduce two auxiliary functions that express the ADMM updates in a form amenable to convergence analysis.
For clarity, we summarize the notation used throughout this section. Let \(S:=NR\) be the total number of axial slices. The indices \(n=1,\dots,S\), \(j=1,\dots,N\), and \(r=1,\dots,R\) refer to individual axial slices, compressed shots, and local slices within a shot, respectively, with \(n=(j-1)R+r\). We stack the slice variables into full vectors \[\boldsymbol{u} := [\boldsymbol{u}_1^\top,\ldots,\boldsymbol{u}_S^\top]^\top \in \mathbb{R}^{pS},\qquad \boldsymbol{v} := [\boldsymbol{v}_1^\top,\ldots,\boldsymbol{v}_S^\top]^\top \in \mathbb{R}^{pS},\qquad \boldsymbol{d} := [\boldsymbol{d}_1^\top,\ldots,\boldsymbol{d}_S^\top]^\top \in \mathbb{R}^{pS},\] and write \(\boldsymbol{u}^k\), \(\boldsymbol{v}^k\), and \(\boldsymbol{d}^k\) for these stacked vectors at the \(k\)th outer iteration. The \(n\)th slice variables at iteration \(k\) are denoted \(\boldsymbol{u}_n^k\), \(\boldsymbol{v}_n^k\), and \(\boldsymbol{d}_n^k\). For each shot \(j\), the corresponding \(R\) slice variables are collected into the shot-wise vector \(\boldsymbol{v}^{(j)}\), and we write \(\boldsymbol{v}^{(j),k}\) for its value at the \(k\)th iteration. In the cyclic GS sweep, each slice variable \(\boldsymbol{u}_n\) is treated as one coordinate block, which we refer to as a slice block; and \(\boldsymbol{u}_{n,-}^k\) and \(\boldsymbol{u}_{n,+}^k\) denote the previous and next neighboring slices used when updating \(\boldsymbol{u}_n\), as defined in 19 . The notation \(\boldsymbol{U}^{k,n}\) denotes the intermediate stacked vector after the first \(n\) slice blocks have been updated in the \(k\)th GS sweep. Both \(\boldsymbol{v}^{(j)}\) and \(\boldsymbol{U}^{k,n}\) are defined explicitly in 25 and 31 , respectively. If \(\boldsymbol{u}^\star\) is an accumulation point, \(\boldsymbol{u}_n^\star\) denotes its \(n\)th slice. All notation is summarized in Table 1.
| Notation | Meaning |
|---|---|
| \(n\) | Global axial-slice index, \(n=1,\dots,S\) |
| \(j\) | Compressed-shot index, \(j=1,\dots,N\) |
| \(r\) | Local slice index within the \(j\)th compressed shot, \(r=1,\dots,R\) |
| \(\bm u_n^k,\bm v_n^k,\bm d_n^k\) | The \(n\)th slice variable at the \(k\)th outer iteration |
| \(\bm u^k,\bm v^k,\bm d^k\) | Stacked vectors at the \(k\)th outer iteration |
| \(\bm v^{(j)}\) | Shot-wise vector collecting the \(R\) slices associated with the \(j\)th shot |
| \(\bm v^{(j),k}\) | Value of the shot-wise vector \(\bm v^{(j)}\) at the \(k\)th iteration |
| \(\bm u_{n,-}^k,\bm u_{n,+}^k\) | Previous and next neighboring slice used when updating the \(n\)th slice |
| \(\bm U^{k,n}\) | Stacked vector after the first \(n\) slices have been updated in the \(k\)th GS sweep |
| \(\bm u_n^\star\) | The \(n\)th slice of an accumulation point \(\bm u^\star\) |
Recall that for each compressed shot \(j=1,\dots,N\), we collect the corresponding \(R\) slice variables into the shot-wise vector \[\label{eq:v94j} \boldsymbol{v}^{(j)} := [\boldsymbol{v}_{(j-1)R+1}^\top,\ldots,\boldsymbol{v}_{jR}^\top]^\top \in\mathbb{R}^{pR},\tag{25}\] and define the shot-wise measurement matrix \(\Phi := [\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_R]\in\mathbb{R}^{p\times pR}.\) Since the same set of masks is applied in every shot, the global forward operator takes the block-diagonal form \[A:=I_N\otimes \Phi\in\mathbb{R}^{pN\times pS},\] where \(\otimes\) denotes the Kronecker product. Stacking all compressed measurements gives \[\boldsymbol{b} := [\boldsymbol{b}_1^\top,\ldots,\boldsymbol{b}_N^\top]^\top \in\mathbb{R}^{pN}.\] The summation in the data term of 15 can then be equivalently written in matrix-vector form as \[\frac{1}{2}\|A\boldsymbol{v}-\boldsymbol{b}\|_2^2 = \frac{1}{2}\sum_{j=1}^{N}\left\|\boldsymbol{b}_j-\sum_{r=1}^{R}\boldsymbol{\phi}_r\boldsymbol{v}_{(j-1)R+r}\right\|_2^2.\]
For the axial coupling term, define the periodic first-order difference operator \[D:\mathbb{R}^{pS}\to\mathbb{R}^{pS},\qquad (D\boldsymbol{u})_n:=\boldsymbol{u}_n-\boldsymbol{u}_{n-1},\qquad \boldsymbol{u}_0:=\boldsymbol{u}_S.\] Then \(\|D\boldsymbol{u}\|_2^2=\sum_{n=1}^{S}\|\boldsymbol{u}_n-\boldsymbol{u}_{n-1}\|_2^2.\) We adopt periodic boundary conditions here for analytical convenience; alternative nonperiodic boundary treatments can be handled similarly.
The slice-based and axial-coupled models can be written in the following unified formulation with \(\gamma\ge 0\): \[\label{eq:convex95unified95stacked} \min_{\boldsymbol{v}\in\mathbb{R}^{pS}}\; F_{\gamma}(\boldsymbol{v}) := \frac{1}{2}\|A\boldsymbol{v}-\boldsymbol{b}\|_2^2 +\lambda\sum_{n=1}^{S}\psi(\boldsymbol{v}_n) +\frac{\gamma}{2}\|D\boldsymbol{v}\|_2^2 .\tag{26}\] When \(\gamma=0\), 26 reduces to the slice-based model; when \(\gamma>0\), it corresponds to the axial-coupled model.
Introduce an auxiliary variable \(\boldsymbol{u}\) and impose the constraint \(\boldsymbol{v}=\boldsymbol{u}\). Then 26 can be written as \[\label{eq:admm95form} \min_{\boldsymbol{v},\boldsymbol{u}}\; f(\boldsymbol{v})+g(\boldsymbol{u}) \quad\text{s.t.}\quad \boldsymbol{v}-\boldsymbol{u}=\boldsymbol{0},\tag{27}\] where \[f(\boldsymbol{v}):=\frac{1}{2}\|A\boldsymbol{v}-\boldsymbol{b}\|_2^2, \qquad g(\boldsymbol{u}):= \lambda\sum_{n=1}^{S}\psi(\boldsymbol{u}_n) +\frac{\gamma}{2}\|D\boldsymbol{u}\|_2^2 .\]
Using the scaled dual variable \(\boldsymbol{d}\), the augmented Lagrangian can be written as \[\label{eq:aug95lag95analysis} \mathcal{L}_\rho(\boldsymbol{u},\boldsymbol{v},\boldsymbol{d}) = f(\boldsymbol{v})+g(\boldsymbol{u}) +\frac{\rho}{2}\|\boldsymbol{v}-\boldsymbol{u}+\boldsymbol{d}\|_2^2 -\frac{\rho}{2}\|\boldsymbol{d}\|_2^2 ,\tag{28}\] or equivalently, \[\label{eq:aug95lag95expanded} \mathcal{L}_\rho(\boldsymbol{u},\boldsymbol{v},\boldsymbol{d}) = f(\boldsymbol{v})+g(\boldsymbol{u}) +\rho\langle \boldsymbol{v}-\boldsymbol{u},\boldsymbol{d}\rangle +\frac{\rho}{2}\|\boldsymbol{v}-\boldsymbol{u}\|_2^2 .\tag{29}\]
Throughout this section, \(\{(\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k)\}_{k\ge0}\) denotes the sequence generated by Algorithm 2.
We now proceed with the convergence analysis. We first review relevant concepts from variational analysis [50], then establish a sufficient decrease property for the augmented Lagrangian in Theorem 1. Building on this and under a standard boundedness assumption, Corollary 1 establishes asymptotic regularity of the iterates, and Theorem 2 shows that every accumulation point of the generated sequence is a stationary point of 26 .
Recall that an extended real-valued function \(h:\mathbb{R}^m\to(-\infty,+\infty]\) is called proper if \(\operatorname{dom} h:=\{x\in\mathbb{R}^m: h(x)<+\infty\}\) is nonempty and \(h(x)>-\infty\) for all \(x\). Additionally, a proper function \(h\) is said to be closed if it is lower semi-continuous. For a proper closed convex function \(h\), its subdifferential at \(x\in\operatorname{dom}h\) is defined by \(\partial h(x) := \{s\in\mathbb{R}^m: h(y)\ge h(x)+\langle s,y-x\rangle,\;\forall y\in\mathbb{R}^m\}.\) When \(h\) is differentiable, we write \(\nabla h\) for its gradient, in which case \(\partial h(x) =\{\nabla h(x)\}.\) A function \(h:\mathbb{R}^m\to(-\infty,+\infty]\) is called \(\kappa\)-weakly convex if \(h(\cdot)+\frac{\kappa}{2}\|\cdot\|_2^2\) is convex for a constant \(\kappa\ge0\), and its subdifferential is defined as \(\partial h(x) := \partial\left(h+\frac{\kappa}{2}\|\cdot\|_2^2\right)(x)-\kappa x.\)
Assumption 1. \(\psi:\mathbb{R}^p\to(-\infty,+\infty]\) is proper, closed, and \(\kappa\)-weakly convex.
The Lipschitz constant of \(\nabla f\) is given by \(L_f:=\|A\|_2^2,\) which equals the largest eigenvalue of \(A^\top A\). With these ingredients in place, the following theorem establishes a sufficient decrease property for the augmented Lagrangian along the iterates of Algorithm 2.
Theorem 1. Suppose Assumption 1 holds and the penalty parameter satisfies \(\rho>2L_f\) and \(\rho+2\gamma>\kappa\lambda\). Then for all \(k\ge 1\), \[\begin{align} &\mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^{k+1},\boldsymbol{d}^{k+1}) - \mathcal{L}_\rho(\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k) \notag\\ &\le -\frac{\rho+2\gamma-\kappa\lambda}{2}\|\boldsymbol{u}^{k+1}-\boldsymbol{u}^k\|_2^2 -\frac{\rho}{2}\|\boldsymbol{v}^{k+1}-\boldsymbol{v}^k\|_2^2 -\left(\frac{1}{2}-\frac{L_f}{\rho}\right) \|A(\boldsymbol{v}^{k+1}-\boldsymbol{v}^k)\|_2^2 . \label{eq:sufficient95decrease} \end{align}\qquad{(1)}\]
Proof. We estimate the changes of the augmented Lagrangian along the three updates.
Fixing \(\boldsymbol{v}^k\) and \(\boldsymbol{d}^k\), we begin with the \(\boldsymbol{u}\)-update. For the \(n\)th slice update \((n=1, \cdots, S)\), we define the single-slice objective \[\begin{align} Q_n^k(\boldsymbol{z}) := \lambda\psi(\boldsymbol{z}) +\frac{\rho}{2}\|\boldsymbol{z}-(\boldsymbol{v}_n^k+\boldsymbol{d}_n^k)\|_2^2 +\frac{\gamma}{2}\|\boldsymbol{z}-\boldsymbol{u}_{n,-}^k\|_2^2 +\frac{\gamma}{2}\|\boldsymbol{z}-\boldsymbol{u}_{n,+}^k\|_2^2 , \label{eq:block95objective} \end{align}\tag{30}\] where \(\boldsymbol{u}_{n,-}^k\) and \(\boldsymbol{u}_{n,+}^k\) are defined in 19 . To handle the GS sweep compactly, we introduce \(\boldsymbol{U}^{k,n}\) to denote the intermediate iterate of the stacked vector after the first \(n\) slices have been updated in the \(k\)th ADMM iteration, that is, \[\label{eq:U} \boldsymbol{U}^{k,n} := [(\boldsymbol{u}_1^{k+1})^\top,\ldots,(\boldsymbol{u}_n^{k+1})^\top, (\boldsymbol{u}_{n+1}^{k})^\top,\ldots,(\boldsymbol{u}_S^{k})^\top]^\top,\tag{31}\] with \(\boldsymbol{U}^{k,0}:=\boldsymbol{u}^k\) and \(\boldsymbol{U}^{k,S}:=\boldsymbol{u}^{k+1}\). At the time of the \(n\)th slice update, \(Q_n^k\) contains exactly the terms in \(\mathcal{L}_\rho(\cdot,\boldsymbol{v}^k,\boldsymbol{d}^k)\) that depend on the slice \(\boldsymbol{u}_n\), with all other slices fixed at their most recently computed values. The remaining terms in \(\mathcal{L}_\rho\) are independent of \(\boldsymbol{z}\) and cancel in the difference, which implies \[\mathcal{L}_\rho(\boldsymbol{U}^{k,n},\boldsymbol{v}^k,\boldsymbol{d}^k) - \mathcal{L}_\rho(\boldsymbol{U}^{k,n-1},\boldsymbol{v}^k,\boldsymbol{d}^k) = Q_n^k(\boldsymbol{u}_n^{k+1})-Q_n^k(\boldsymbol{u}_n^k).\] Under Assumption 1 and the the condition that \(\rho+2\gamma>\kappa\lambda\), \(Q_n^k\) is \((\rho+2\gamma-\kappa\lambda)\)-strongly convex. Since the \(n\)th slice update minimizes \(Q_n^k\), we have \[Q_n^k(\boldsymbol{u}_n^{k+1})-Q_n^k(\boldsymbol{u}_n^k) \le -\frac{\rho+2\gamma-\kappa\lambda}{2}\|\boldsymbol{u}_n^{k+1}-\boldsymbol{u}_n^k\|_2^2 .\] Consequently, \[\mathcal{L}_\rho(\boldsymbol{U}^{k,n},\boldsymbol{v}^k,\boldsymbol{d}^k) - \mathcal{L}_\rho(\boldsymbol{U}^{k,n-1},\boldsymbol{v}^k,\boldsymbol{d}^k) \le -\frac{\rho+2\gamma-\kappa\lambda}{2}\|\boldsymbol{u}_n^{k+1}-\boldsymbol{u}_n^k\|_2^2 .\] Summing over \(n=1,\dots,S\) yields \[\label{eq:u95decrease} \mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^k,\boldsymbol{d}^k) - \mathcal{L}_\rho(\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k) \le -\frac{\rho+2\gamma-\kappa\lambda}{2}\|\boldsymbol{u}^{k+1}-\boldsymbol{u}^k\|_2^2 .\tag{32}\]
Next, we estimate the decrease from the \(\boldsymbol{v}\)-update. Let \[\Delta\boldsymbol{v}^{k+1}:=\boldsymbol{v}^{k+1}-\boldsymbol{v}^k .\] For fixed \(\boldsymbol{u}^{k+1}\) and \(\boldsymbol{d}^k\), define the \(\boldsymbol{v}\)-dependent quadratic function \[q_k(\boldsymbol{v}) := \frac{1}{2}\|A\boldsymbol{v}-\boldsymbol{b}\|_2^2 + \rho\langle \boldsymbol{v}-\boldsymbol{u}^{k+1},\boldsymbol{d}^k\rangle + \frac{\rho}{2}\|\boldsymbol{v}-\boldsymbol{u}^{k+1}\|_2^2 .\] By the expanded form 29 , the difference, \[\mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^{k+1},\boldsymbol{d}^k) - \mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^k,\boldsymbol{d}^k),\] is equal to \(q_k(\boldsymbol{v}^{k+1})-q_k(\boldsymbol{v}^k)\). Since \(q_k\) is quadratic with Hessian \(A^\top A+\rho I\), we have \[\begin{align} q_k(\boldsymbol{v}^{k+1})-q_k(\boldsymbol{v}^k) =&\; \left\langle \nabla q_k(\boldsymbol{v}^{k+1}),\Delta\boldsymbol{v}^{k+1}\right\rangle -\frac{1}{2}\|A\Delta\boldsymbol{v}^{k+1}\|_2^2 -\frac{\rho}{2}\|\Delta\boldsymbol{v}^{k+1}\|_2^2 . \label{eq:v95step95expand} \end{align}\tag{33}\] The optimality condition of the \(\boldsymbol{v}\)-update is \[\nabla q_k(\boldsymbol{v}^{k+1}) = A^\top(A\boldsymbol{v}^{k+1}-\boldsymbol{b}) + \rho(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}+\boldsymbol{d}^k) = \boldsymbol{0}.\] Therefore, \[\label{eq:v95decrease} \mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^{k+1},\boldsymbol{d}^k) - \mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^k,\boldsymbol{d}^k) = -\frac{1}{2}\|A\Delta\boldsymbol{v}^{k+1}\|_2^2 -\frac{\rho}{2}\|\Delta\boldsymbol{v}^{k+1}\|_2^2 .\tag{34}\]
It remains to control the change caused by the dual update. From the \(\boldsymbol{v}\)-optimality condition and \[\boldsymbol{d}^{k+1}=\boldsymbol{d}^k+\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1},\] we obtain \[\label{eq:dual95relation} A^\top(A\boldsymbol{v}^{k+1}-\boldsymbol{b})+\rho\boldsymbol{d}^{k+1}=\boldsymbol{0},\tag{35}\] or equivalently, \[\label{eq:dual95relation2} \boldsymbol{d}^{k+1} = -\frac{1}{\rho} A^\top(A\boldsymbol{v}^{k+1}-\boldsymbol{b}).\tag{36}\] Applying this relation at two consecutive iterations gives, for \(k\ge 1\), \[\label{eq:difference95dual} \boldsymbol{d}^{k+1}-\boldsymbol{d}^k = -\frac{1}{\rho} A^\top A(\boldsymbol{v}^{k+1}-\boldsymbol{v}^k).\tag{37}\] Using the expanded form 29 , we also have \[\mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^{k+1},\boldsymbol{d}^{k+1}) - \mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^{k+1},\boldsymbol{d}^k) = \rho\left\langle \boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}, \boldsymbol{d}^{k+1}-\boldsymbol{d}^k \right\rangle .\] It follows from \(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}=\boldsymbol{d}^{k+1}-\boldsymbol{d}^k\) together with 37 that \[\begin{align} \notag &\mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^{k+1},\boldsymbol{d}^{k+1}) - \mathcal{L}_\rho(\boldsymbol{u}^{k+1},\boldsymbol{v}^{k+1},\boldsymbol{d}^k)\\ =& \rho\|\boldsymbol{d}^{k+1}-\boldsymbol{d}^k\|_2^2 = \frac{1}{\rho} \|A^\top A(\boldsymbol{v}^{k+1}-\boldsymbol{v}^k)\|_2^2 \le \frac{L_f}{\rho} \|A(\boldsymbol{v}^{k+1}-\boldsymbol{v}^k)\|_2^2 . \label{eq:dual95increase} \end{align}\tag{38}\]
Combining the estimates for the \(\boldsymbol{u}\)-step in 32 , the \(\boldsymbol{v}\)-step in 34 , and the dual update in 38 gives ?? . ◻
Assumption 2. The function \(\psi\) is bounded below.
This assumption holds for commonly used nonnegative regularizers, such as Tikhonov and TV regularization.
Lemma 1. Suppose Assumption 2 holds and \(\rho>2L_f\). Then the augmented Lagrangian sequence \(\{\mathcal{L}_\rho(\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k)\}_{k\ge 1}\) is bounded below.
Proof. Using the relation at the \(k\)th iteration, we obtain \[\begin{align} \notag \mathcal{L}_\rho(\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k) =& f(\boldsymbol{v}^k)+g(\boldsymbol{u}^k) +\frac{\rho}{2}\|\boldsymbol{v}^k-\boldsymbol{u}^k+\boldsymbol{d}^k\|_2^2 -\frac{\rho}{2}\|\boldsymbol{d}^k\|_2^2 \\ \ge& f(\boldsymbol{v}^k)+g(\boldsymbol{u}^k)-\frac{1}{2\rho}\|A^\top(A\boldsymbol{v}^{k}-\boldsymbol{b})\|_2^2.\label{ineq:lower95bound} \end{align}\tag{39}\] Since \(f(\boldsymbol{v})=\frac{1}{2}\|A\boldsymbol{v}-\boldsymbol{b}\|_2^2\), it is straightforward that \[\|\nabla f(\boldsymbol{v})\|_2^2 = \|A^\top(A\boldsymbol{v}-\boldsymbol{b})\|_2^2 \le \|A\|_2^2\|A\boldsymbol{v}-\boldsymbol{b}\|_2^2 = 2L_f f(\boldsymbol{v}).\] We further plug the expression of \(g(\boldsymbol{u}^k)\) into 39 and obtain \[\mathcal{L}_\rho(\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k) \ge \left(1-\frac{L_f}{\rho}\right)f(\boldsymbol{v}^k)+g(\boldsymbol{u}^k) =\left(1-\frac{L_f}{\rho}\right)f(\boldsymbol{v}^k)+\lambda\sum_{n=1}^{S}\psi(\boldsymbol{u}_n^k) +\frac{\gamma}{2}\|D\boldsymbol{u}^k\|_2^2.\] Since \(\rho>2L_f\), we have \(1-L_f/\rho>0\). The boundedness of \(\psi\) in Assumption 2 guarantees that \(\{\mathcal{L}_\rho(\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k)\}_{k\ge1}\) is bounded below. ◻
Corollary 1. Suppose Assumptions 1–2 hold, \(\rho>2L_f,\) and \(\rho+2\gamma>\kappa\lambda\). Then \(\sum_{k=1}^{\infty}\|\boldsymbol{u}^{k+1}-\boldsymbol{u}^k\|_2^2<\infty\) and \(\sum_{k=1}^{\infty}\|\boldsymbol{v}^{k+1}-\boldsymbol{v}^k\|_2^2<\infty\). Consequently, \[\|\boldsymbol{u}^{k+1}-\boldsymbol{u}^k\|_2\to 0,\qquad \|\boldsymbol{v}^{k+1}-\boldsymbol{v}^k\|_2\to 0,\qquad \|\boldsymbol{d}^{k+1}-\boldsymbol{d}^k\|_2\to 0,\] and the primal residual satisfies \(\|\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}\|_2\to 0\).
Proof. Since \(\rho>2L_f\) and \(\rho+2\gamma>\kappa\lambda\), the coefficients \(\frac{1}{2}-\frac{L_f}{\rho}\) and \(\frac{\rho+2\gamma-\kappa\lambda}{2}\) in ?? are positive. Summing ?? over \(k\) and using Lemma 1 gives \(\sum_{k=1}^{\infty}\|\boldsymbol{u}^{k+1}-\boldsymbol{u}^k\|_2^2<\infty\) and \(\sum_{k=1}^{\infty}\|\boldsymbol{v}^{k+1}-\boldsymbol{v}^k\|_2^2<\infty\). Therefore, \(\|\boldsymbol{u}^{k+1}-\boldsymbol{u}^k\|_2\to 0\) and \(\|\boldsymbol{v}^{k+1}-\boldsymbol{v}^k\|_2\to 0\). The relation 37 implies \(\|\boldsymbol{d}^{k+1}-\boldsymbol{d}^k\|_2\to 0\). Finally, by the dual update, \(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}=\boldsymbol{d}^{k+1}-\boldsymbol{d}^k\), and hence the primal residual also converges to zero. ◻
Lemma 2. Suppose Assumption 1 holds. For every \(k\ge 1\), there exist vectors \(\boldsymbol{e}_u^{k+1}\), \(\boldsymbol{e}_v^{k+1}\), and \(\boldsymbol{e}_d^{k+1}\) such that \[\begin{align} \boldsymbol{e}_u^{k+1} &\in \partial g(\boldsymbol{u}^{k+1}) -\rho(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}+\boldsymbol{d}^{k+1}), \label{eq:relative95error95u}\\ \boldsymbol{e}_v^{k+1} &= \nabla f(\boldsymbol{v}^{k+1}) +\rho(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}+\boldsymbol{d}^{k+1}), \label{eq:relative95error95v}\\ \boldsymbol{e}_d^{k+1} &= \rho(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}), \label{eq:relative95error95d} \end{align}\] {#eq: sublabel=eq:eq:relative95error95u,eq:eq:relative95error95v,eq:eq:relative95error95d} and \[\label{eq:relative95error95bound} \|\boldsymbol{e}_u^{k+1}\|_2+\|\boldsymbol{e}_v^{k+1}\|_2+\|\boldsymbol{e}_d^{k+1}\|_2 \le C_{\rm rel} \left( \|\boldsymbol{u}^{k+1}-\boldsymbol{u}^k\|_2 + \|\boldsymbol{v}^{k+1}-\boldsymbol{v}^k\|_2 \right),\qquad{(2)}\] where \(C_{\rm rel}>0\) is independent of \(k\). Consequently, under the assumptions of Corollary 1, \[\|\boldsymbol{e}_u^{k+1}\|_2+\|\boldsymbol{e}_v^{k+1}\|_2+\|\boldsymbol{e}_d^{k+1}\|_2\to 0 .\]
Proof. Let \(\Delta\boldsymbol{u}^{k+1}:=\boldsymbol{u}^{k+1}-\boldsymbol{u}^k\), \(\Delta\boldsymbol{v}^{k+1}:=\boldsymbol{v}^{k+1}-\boldsymbol{v}^k\), \(\Delta\boldsymbol{d}^{k+1}:=\boldsymbol{d}^{k+1}-\boldsymbol{d}^k\). The optimality condition of the \(\boldsymbol{v}\)-update gives \(\nabla f(\boldsymbol{v}^{k+1}) + \rho(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}+\boldsymbol{d}^k) = \boldsymbol{0} .\) Thus, with \(\boldsymbol{e}_v^{k+1}\) defined by ?? , \[\label{eq:ev95bound} \|\boldsymbol{e}_v^{k+1}\|_2 = \rho\|\Delta\boldsymbol{d}^{k+1}\|_2 .\tag{40}\] Moreover, the dual update gives \(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}=\Delta\boldsymbol{d}^{k+1}.\) Thus, with \(\boldsymbol{e}_d^{k+1}\) defined by ?? , \[\label{eq:ed95bound} \|\boldsymbol{e}_d^{k+1}\|_2 = \rho\|\Delta\boldsymbol{d}^{k+1}\|_2 .\tag{41}\]
It remains to estimate the \(\boldsymbol{u}\)-component. The optimality condition of 20 implies that for every slice \(n\), there exists \(\boldsymbol{s}_n^{k+1}\in\partial\psi(\boldsymbol{u}_n^{k+1})\) such that \[\label{eq:u95block95opt95relative} \lambda\boldsymbol{s}_n^{k+1} + \rho(\boldsymbol{u}_n^{k+1}-\boldsymbol{v}_n^k-\boldsymbol{d}_n^k) + \gamma(\boldsymbol{u}_n^{k+1}-\boldsymbol{u}_{n,-}^{k}) + \gamma(\boldsymbol{u}_n^{k+1}-\boldsymbol{u}_{n,+}^{k}) = \boldsymbol{0} .\tag{42}\] Define \[\boldsymbol{q}_n^{k+1} := \lambda\boldsymbol{s}_n^{k+1} + \gamma(2\boldsymbol{u}_n^{k+1}-\boldsymbol{u}_{n-1}^{k+1}-\boldsymbol{u}_{n+1}^{k+1}), \qquad n=1,\dots,S,\] where \(\boldsymbol{u}_0^{k+1}=\boldsymbol{u}_S^{k+1}\) and \(\boldsymbol{u}_{S+1}^{k+1}=\boldsymbol{u}_1^{k+1}\), and let \(\boldsymbol{q}^{k+1}:= [(\boldsymbol{q}_1^{k+1})^\top,\ldots,(\boldsymbol{q}_S^{k+1})^\top]^\top .\) By the definition of \(g\), we have \(\boldsymbol{q}^{k+1}\in\partial g(\boldsymbol{u}^{k+1}).\) Now define \[\label{eq:eu95def95relative} \boldsymbol{e}_u^{k+1} := \boldsymbol{q}^{k+1} - \rho(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}+\boldsymbol{d}^{k+1}).\tag{43}\] Then ?? holds. Moreover, by 42 , the \(n\)th slice of \(\boldsymbol{e}_u^{k+1}\) satisfies \[\label{eq:eu95identity95relative} \boldsymbol{e}_{u,n}^{k+1} = -\rho(\boldsymbol{v}_n^{k+1}-\boldsymbol{v}_n^k) -\rho(\boldsymbol{d}_n^{k+1}-\boldsymbol{d}_n^k) +\gamma(\boldsymbol{u}_{n,-}^{k}-\boldsymbol{u}_{n-1}^{k+1}) +\gamma(\boldsymbol{u}_{n,+}^{k}-\boldsymbol{u}_{n+1}^{k+1}).\tag{44}\] The last two terms come from the mixed old and new neighboring slices in one cyclic Gauss–Seidel sweep. By 19 , \[\left( \sum_{n=1}^{S} \left\| (\boldsymbol{u}_{n,-}^{k}-\boldsymbol{u}_{n-1}^{k+1}) + (\boldsymbol{u}_{n,+}^{k}-\boldsymbol{u}_{n+1}^{k+1}) \right\|_2^2 \right)^{1/2} \le 2\|\Delta\boldsymbol{u}^{k+1}\|_2 .\] Therefore, \[\label{eq:eu95bound} \|\boldsymbol{e}_u^{k+1}\|_2 \le 2\gamma\|\Delta\boldsymbol{u}^{k+1}\|_2 + \rho\|\Delta\boldsymbol{v}^{k+1}\|_2 + \rho\|\Delta\boldsymbol{d}^{k+1}\|_2 .\tag{45}\]
Finally, by 37 , \(\Delta\boldsymbol{d}^{k+1} = -\frac{1}{\rho} A^\top A\Delta\boldsymbol{v}^{k+1},\) and hence \(\rho\|\Delta\boldsymbol{d}^{k+1}\|_2 \le L_f\|\Delta\boldsymbol{v}^{k+1}\|_2 .\) Combining this estimate with 40 , 41 , and 45 gives ?? with any \(C_{\rm rel}\ge 2\gamma+\rho+3L_f .\) The final claim follows from Corollary 1. ◻
Theorem 2. Suppose Assumptions 1–2 hold, \(\rho>2L_f\), and \(\rho+2\gamma>\kappa\lambda\). If the sequence \(\{(\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k)\}\) generated by Algorithm 2 is bounded, then every accumulation point \((\boldsymbol{u}^\star,\boldsymbol{v}^\star,\boldsymbol{d}^\star)\) satisfies \(\boldsymbol{v}^\star=\boldsymbol{u}^\star\) and \[\label{eq:kkt95limit} \boldsymbol{0}\in \nabla f(\boldsymbol{v}^\star)+\partial g(\boldsymbol{v}^\star).\qquad{(3)}\] Consequently, \(\boldsymbol{v}^\star\) is a stationary point of 26 .
Proof. Let \((\boldsymbol{u}^\star,\boldsymbol{v}^\star,\boldsymbol{d}^\star)\) be an accumulation point of \(\{(\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k)\}\), which exists by boundedness and the Bolzano-Weierstrass theorem. Take a subsequence converging to it; without loss of generality, we index this subsequence by \(k\) for simplicity. By Corollary 1, the shifted subsequence \((\boldsymbol{u}^{k+1},\boldsymbol{v}^{k+1},\boldsymbol{d}^{k+1})\) converges to the same limit. Since \(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}\to \boldsymbol{0}\), we have \(\boldsymbol{v}^\star=\boldsymbol{u}^\star\).
The optimality condition of the \(\boldsymbol{v}\)-update gives \(\nabla f(\boldsymbol{v}^{k+1})+\rho\boldsymbol{d}^{k+1}=\boldsymbol{0}.\) Passing to the limit yields \[\label{eq:v95limit95opt} \nabla f(\boldsymbol{v}^\star)+\rho\boldsymbol{d}^\star=\boldsymbol{0}.\tag{46}\] By Lemma 2, there exist vectors \(\boldsymbol{e}_u^{k+1}\), \(\boldsymbol{e}_v^{k+1}\), and \(\boldsymbol{e}_d^{k+1}\) satisfying ?? –?? , such that \(\|\boldsymbol{e}_u^{k+1}\|_2+\|\boldsymbol{e}_v^{k+1}\|_2+\|\boldsymbol{e}_d^{k+1}\|_2\to 0.\) From ?? , there exists \(\boldsymbol{q}^{k+1}\in\partial g(\boldsymbol{u}^{k+1})\) such that \(\boldsymbol{e}_u^{k+1} = \boldsymbol{q}^{k+1} - \rho(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}+\boldsymbol{d}^{k+1}).\) Since \(\boldsymbol{e}_u^{k+1}\to\boldsymbol{0}\), \(\boldsymbol{v}^{k+1}-\boldsymbol{u}^{k+1}\to\boldsymbol{0}\), and \(\boldsymbol{d}^{k+1}\to\boldsymbol{d}^\star\), we have \(\boldsymbol{q}^{k+1}\to \rho\boldsymbol{d}^\star .\) Under Assumption 1, the function \(g\) is proper, closed, and weakly convex, so its subdifferential is closed under limits. Since \(\boldsymbol{q}^{k+1}\in\partial g(\boldsymbol{u}^{k+1})\), \(\boldsymbol{u}^{k+1}\to\boldsymbol{u}^\star\), and \(\boldsymbol{q}^{k+1}\to\rho\boldsymbol{d}^\star\), we obtain \[\label{eq:u95limit95opt} \rho\boldsymbol{d}^\star\in\partial g(\boldsymbol{u}^\star).\tag{47}\]
Since \(\boldsymbol{u}^\star=\boldsymbol{v}^\star\), combining 46 and 47 gives \[\boldsymbol{0}\in \nabla f(\boldsymbol{v}^\star)+\partial g(\boldsymbol{v}^\star).\] Therefore, \(\boldsymbol{v}^\star\) is a stationary point of 26 . ◻
The boundedness of the generated sequence \(\{\boldsymbol{u}^k,\boldsymbol{v}^k,\boldsymbol{d}^k\}\) is assumed rather than proved. A sufficient condition that guarantees boundedness is coercivity [51] of \(\psi,\) which ensures that the augmented Lagrangian sublevel sets are compact. Both Tikhonov and TV regularizations satisfy this condition. Since coercivity may be difficult to verify for a general function \(\psi,\) we state boundedness directly as a hypothesis, which can be monitored empirically during the iteration.
The weakly convex assumption on \(\psi\) provides a natural middle ground between convex and fully nonconvex regularization. Classical convergence guarantees [31], [47] require convexity, which can be restrictive in practice, as some useful regularization models in image reconstruction are nonconvex. Weak convexity retains enough structure to support the convergence analysis above while accommodating a broader class of regularizers. This connects to recent efforts to extend PnP convergence theory to nonconvex settings: for example, Hurault et al. [52] propose a proximal denoising step that, in certain cases, coincides with the exact proximal operator of a possibly nonconvex regularizer, and Shoushtari et al. [53] establish nonconvex convergence guarantees for PnP-ADMM under distribution shift, focusing on the case where a learned denoiser is applied outside its training distribution.
The convergence analysis in this section is variational in nature: it assumes that the denoising step is induced by a regularization functional \(\psi\). Another line of PnP research establishes convergence by imposing structural assumptions directly on the denoiser. For example, convergence results for PnP methods with black-box denoisers may require denoiser-based assumptions, such as nonexpansiveness, contractiveness, or averagedness of the denoiser itself or of a related residual operator [48], [49]. Such results often lead to fixed-point convergence guarantees rather than convergence for a variational model. In contrast, our analysis assumes that \(\psi\) is weakly convex and proves subsequential convergence for the corresponding variational formulation, and the result applies directly to Tikhonov and TV. For off-the-shelf denoisers such as BM3D, DnCNN, FFDNet, or DRUNet, the iteration is generally no longer equivalent to minimizing a fixed objective. Establishing or verifying such denoiser-based assumptions within the present framework is beyond the scope of this work; for these denoisers, we report empirical behavior only. Developing a rigorous convergence theory for black-box denoisers within the proposed axial-coupled PnP scheme remains an open problem for future work.
4 Experiments↩︎
In this section, we evaluate the proposed reconstruction framework on synthetic and real zebrafish-heart CS-LSM data. The synthetic experiment provides a controlled setting with a known reference volume, whereas the real-data experiment evaluates the method under realistic acquisition conditions using a physical imaging system. The experiments evaluate the effects of denoiser choice, axial coupling, and the compression ratio. All experiments were implemented in Python and run on an Azure Machine Learning Standard_E4ds_v4 CPU instance with 4 cores and 32 GB RAM.
4.1 Experimental setup↩︎
4.1.0.1 Datasets
For the synthetic experiment, we use a zebrafish-heart reference volume derived from LSM data acquired by our in-house LSM system [54]. The underlying dataset was obtained using retrospective synchronization, which aligns image sequences from different axial slices and cardiac cycles. The selected reference volume (\(200\times200\times120\)-voxel) corresponds to a representative time point, containing both the atrium and ventricle with approximately \(300\) nuclei near the chamber surfaces. Compressed measurements are generated from this reference using the forward model in 1 with known random binary masks, consistent with the DMD coding strategy in CS-LSM. The compression ratio is set to \(R=5\) unless otherwise stated. Results for varying compression ratios are reported in Section 4.4. No additional measurement noise is added in the synthetic experiments, so that the comparison isolates the effects of the inverse solver, the denoising prior, and the axial coupling term.
For the real-data experiment, we use CS-LSM measurements and the corresponding sensing masks acquired by the DMD-based light-sheet platform described in [29]. The real dataset consists of \(N=10\) compressed camera shots, each of size \(404\times404\) pixels, acquired with compression ratio \(R=15\), thus the reconstructed volume contains \(S=150\) axial slices. These measurements are collected under imaging conditions comparable to those of the synthetic data in zebrafish-heart, but a corresponding volumetric ground truth is not available. To focus on reconstruction performance, we compare different models and denoisers directly on the acquired measurements without additional pre- or post-processing. While addressing effects such as rolling-shutter distortion may further improve image quality, we refer readers to our related work [29] for a detailed discussion of these issues.
4.1.0.2 Evaluation metrics
To assess reconstruction performance on synthetic data, we adopt two standard image-quality metrics: Peak Signal-to-Noise Ratio (PSNR) [55] and Structural Similarity Index Measure (SSIM) [56]. PSNR measures pixel-wise reconstruction fidelity, with larger values indicating smaller error. SSIM measures structural similarity between each reconstructed slice and the corresponding reference, following [56], with values closer to 1 indicating greater similarity. Both metrics are computed per slice and averaged over all valid slices. For real data, we report visual comparisons only.
4.1.0.3 Parameter tuning
The PnP-ADMM framework has three hyperparameters: the ADMM penalty parameter \(\rho\), the parameter \(\lambda\) controlling the strength of the regularization or denoising step, and the axial coupling weight \(\gamma\), with \(\gamma=0\) for the slice-based model. For each denoiser and reconstruction model, we tune these parameters on the synthetic dataset using Bayesian optimization [57], [58], with negative slice-averaged PSNR as the objective. Minimizing this objective is equivalent to maximizing the average reconstruction quality over the axial stack. Each Bayesian-optimization run is limited to \(50\) objective evaluations over predefined search ranges, and the optimized parameters are then used to rerun the corresponding reconstruction. The optimized hyperparameters for all denoiser and model combinations are reported in Table 2. For the real data, no separate parameter search is performed, as an exactly matched volumetric ground truth is unavailable; instead, parameters selected on the synthetic dataset are reused, given similar imaging conditions.
The optimal parameter \(\rho\) obtained via Bayesian optimization is notably smaller than the value required by Theorem 2. This is not a contradiction, as Theorem 2 provides a sufficient (not necessary) condition for convergence of the axial-coupled PnP-ADMM scheme, and the true convergence region may be considerably larger. As shown in Table 2, smaller values of \(\rho\) work well empirically and tend to yield better reconstruction quality, suggesting that the theoretical bound is conservative. Such a gap between sufficient conditions and practical performance is common in convergence analysis.
4.1.0.4 Stopping criteria
For synthetic data, the PnP-ADMM iterations are terminated when the relative change in iterates, i.e., \(\|\boldsymbol{v}^{k}-\boldsymbol{v}^{k-1}\|_2/\|\boldsymbol{v}^{k-1}\|_2\), falls below \(10^{-3}\), or when the number of iterations reaches \(100\). For real data, a looser tolerance of \(10^{-2}\) is used to account for measurement noise.
| Denoiser | Slice-based | Axial-coupled | |||
| \(\rho\) | \(\lambda\) | \(\rho\) | \(\lambda\) | \(\gamma\) | |
| Tikhonov | 0.0501 | 0.0198 | 0.0412 | 0.0248 | 0.0097 |
| TV | 0.4999 | 0.0673 | 0.0068 | 0.0018 | 0.0010 |
| BM3D | 0.7828 | 0.0075 | 0.0083 | 0.0001 | 0.0011 |
| DnCNN | 0.0812 | 0.0008 | 0.0156 | 0.0002 | 0.0031 |
| FFDNet | 0.1439 | 0.0014 | 0.0403 | 0.0004 | 0.0010 |
| DRUNet | 0.1921 | 0.0018 | 0.0236 | 0.0002 | 0.0010 |
5pt
4.2 Synthetic data↩︎
4.2.0.1 Quantitative results
| Denoiser | PSNR (dB) \(\uparrow\) | SSIM \(\uparrow\) | Runtime (s) \(\downarrow\) | |||
| Slice | Axial | Slice | Axial | Slice | Axial | |
| Tikhonov | 29.7675 | 33.9182 | 0.7230 | 0.9283 | 0.579 | 1.267 |
| TV | 36.4616 | 40.1035 | 0.9551 | 0.9770 | 1.861 | 3.674 |
| BM3D | 37.8380 | 40.2827 | 0.9617 | 0.9718 | 392.380 | 501.994 |
| DnCNN | 42.0069 | 42.8288 | 0.9751 | 0.9735 | 121.921 | 142.108 |
| FFDNet | 43.1213 | 43.3507 | 0.9859 | 0.9866 | 61.005 | 68.156 |
| DRUNet | 43.2100 | 43.7359 | 0.9866 | 0.9877 | 406.829 | 504.127 |
2.5pt
Table 3 summarizes the results on synthetic data for the slice-based and axial-coupled models. Compared with slice-based reconstructions, axial coupling improves PSNR across all considered denoisers. SSIM likewise improves in most cases, with the only exception being DnCNN, where the difference is negligible. Notably, the improvement is more pronounced for classical priors such as Tikhonov and TV, where slice-based reconstructions are less accurate, but smaller for deep-learning-based denoisers such as FFDNet and DRUNet, whose slice-based reconstructions are already close to the reference.
The gains from axial coupling reflect the multiplexed structure of the CS-LSM inverse problem. Under compressed acquisition, each camera measurement is formed by the masked superposition of \(R\) axial slices. Although the data-fidelity term constrains these slices jointly, the slice-based model does not use neighboring slices as an additional prior. As a result, weak cellular or nuclear signals can be less stably localized along the axial direction and may appear fragmented or inconsistent across neighboring slices. The axial term therefore acts as a weak consistency constraint along \(z\), encouraging neighboring slices to share structural information while avoiding strong axial smoothing. Except for Tikhonov, the selected \(\gamma\) values in Table 2 are on the order of \(10^{-3}\), suggesting that weak axial coupling is effective for most denoisers.
The results reveal a clear dependence on the choice of denoiser. Tikhonov is the fastest method but yields the lowest PSNR and SSIM, while TV provides a stronger classical baseline with only a modest increase in runtime. BM3D improves upon TV in the slice-based setting, although its runtime is substantially higher due to the nonlocal patch search. Among the DL-based denoisers, FFDNet and DRUNet achieve the best reconstruction quality. DRUNet attains the highest PSNR/SSIM in both the slice-based and axial-coupled models, whereas FFDNet achieves comparable accuracy with a much lower computational cost, offering a better overall balance between reconstruction quality and efficiency.
Figure 4: Intensity line profiles from the red zoomed region in Fig. 3. (a) Profile location. (b) Slice-based reconstructions. (c) Axial-coupled reconstructions. The profiles show peak separation and local contrast relative to the ground truth..
4.2.0.2 Qualitative results
In all qualitative visualizations, the ventricle and atrium of zebrafish-heart are labeled as V and A respectively for anatomical reference, and all scale bars represent 30 \(\mu m\). Figure 3 presents visual reconstruction results on the synthetic dataset. For each method, we display slice 15 as a representative axial slice, together with the maximum intensity projection (MIP) along the axial direction. MIP is a standard volume-rendering technique that projects the maximum intensity along viewing ray [59], making it well suited for visualizing bright fluorescent structures in 3D. Here, the slice view assesses local image fidelity at a fixed depth, while the MIP provides a global visual summary of the reconstructed volume, emphasizing volumetric continuity, background suppression, and the recovery of bright cellular or nuclear signals in 3D. We further highlight two zoomed-in regions (red and green) to illustrate key reconstruction challenges: avoiding the merging or distortion of adjacent bright structures and preserving weak signals in low-intensity regions. These features are critical for downstream tasks, including cell detection, segmentation, and 3D/4D tracking in zebrafish-heart imaging. Undetected weak signals, merged neighboring structures, and spurious background artifacts can bias the quantitative analyses of cardiac morphology and contractile motion [24], [54].
The main visual differences across denoisers lie in fine-structure recovery and background cleanliness. Tikhonov recovers the coarse cardiac shape but leaves a noisy background in the slice-based reconstruction and suppresses many weak cellular or nuclear signals. TV gives a cleaner and sharper baseline, but it still oversmooths weak signals and thin structures. BM3D preserves more local detail than TV, although residual low-intensity artifacts remain visible. The DL-based denoisers recover brighter and better-separated cellular or nuclear signals, with FFDNet and DRUNet appearing closest to the ground truth. This is especially visible in the red zoomed region, where FFDNet and DRUNet better preserve the separation and morphology of neighboring bright cellular or nuclear signals. To further quantify this local comparison, Fig. 4 shows intensity profiles as a function of spatial position along a representative oblique line in the red zoomed region. The profiles are plotted against physical distance in \(\mu\)m and globally normalized using the intensity range across all compared images. FFDNet and DRUNet more closely match the ground-truth peak locations and valley structure, whereas Tikhonov and TV tend to suppress peak amplitudes and reduce local contrast.
The effect of axial coupling is more evident in the MIP views, where slice-to-slice inconsistencies accumulate across the axial stack. Compared with the slice-based reconstructions, the axial-coupled results generally show cleaner backgrounds and cardiac boundaries, and more stable recovery of cellular structures across \(z\). The improvement of axial coupling appears in different local regions for different denoisers. For example, BM3D shows improved separation and morphology of adjacent bright signals in the red zoomed region, while FFDNet recovers weak signals in the green zoomed region closer to the ground truth. Comparing Figures 4 (b) and 4 (c), the axial-coupled profiles are generally closer to the ground truth for the stronger denoisers, especially in preserving the two-peak structure and the intervening valley.
Overall, the qualitative results agree with the quantitative comparison in Table 3. FFDNet and DRUNet provide the strongest visual reconstructions, while axial coupling improves the consistency of the reconstructed volume, especially for weak cellular or nuclear signals and for classical priors such as Tikhonov and TV. Together, these results highlight that the choice of denoiser determines the baseline reconstruction quality, while axial coupling provides additional benefit by leveraging inter-slice correlations, particularly when the slice-based prior is relatively weak.
4.3 Real data↩︎
Figure 5 shows real-data reconstruction results using MIPs. Because the real-data evaluation is qualitative, we focus on structural continuity, local contrast, background suppression, and the visibility of cellular or nuclear signals rather than voxel-wise fidelity. The comparison examines practical trade-offs among different denoisers, specifically in terms of weak-signal preservation, background suppression, and the recovery of plausible cardiac morphology.
Among the classical priors, Tikhonov and TV retain the main cardiac region and many diffuse low-intensity structures, but at the cost of noticeable background haze and relatively low local contrast. BM3D reduces background variation and yields a cleaner appearance, though some weak structures may also be attenuated in the process. The DL-based denoisers produce sharper bright signals and higher local contrast, but they tend to emphasize the most prominent structures while suppressing some diffuse weak signals. Among them, FFDNet and DRUNet provide the clearest separation between bright nuclear signals and the background.
Axial coupling primarily improves the real-data MIPs by reducing background haze and enhancing the visibility of bright structures. This effect is most pronounced for the DL-based denoisers, where the axial-coupled reconstructions exhibit sharper and more distinctly separated cellular or nuclear signals compared with slice-based results. However, without a matched ground truth, sharper punctate signals alone do not establish improved reconstruction fidelity; they may reflect genuine contrast improvement, denoiser-induced enhancement, or both. Their practical relevance therefore depends on the downstream task.
The real-data experiment presents additional challenges: the effective sensing operator may deviate from the ideal forward model due to mask misalignment, nonideal modulation, or calibration errors, and the noise statistics may not exactly match the assumed model. Without ground truth, the comparison instead reveals practical visual trade-offs: classical priors tend to retain diffuse low-intensity content but leave stronger background haze, while DL-based denoisers improve local contrast and background suppression at the cost of reduced visibility of weaker diffuse signals. This trade-off highlights a practical advantage of the PnP formulation, namely that the reconstruction prior can be selected according to the needs of the downstream task.
4.4 Effect of compression ratio↩︎
We investigate the dependence of reconstruction quality on the compression ratio \(R\), defined as the number of axial slices multiplexed into one camera exposure. Larger values of \(R\) improve acquisition efficiency by reducing the number of measurements per volume, but also increase the degree of ill-posedness of the inverse problem.
Figure 6 reports PSNR and SSIM on the synthetic dataset for TV and FFDNet, representing classical and DL-based denoisers, respectively. For both methods, reconstruction quality generally decreases as \(R\) increases because each reconstructed volume is supported by fewer measurements. For TV, the axial-coupled model consistently outperforms the slice-based model across the tested compression ratios, especially in SSIM. This indicates that inter-slice correlation provides useful structural support when the slice-wise prior is relatively simple. For FFDNet, the slice-based reconstruction is already strong at low compression ratios, so the additional gain from axial coupling is smaller and not uniformly positive across all \(R\). However, at higher compression ratios like \(R=10\), the axial-coupled FFDNet reconstruction improves both PSNR and SSIM over slice-based result.
Figure 6: Effect of compression ratio \(R\) on synthetic reconstruction quality for TV and FFDNet. Each curve reports slice-averaged PSNR or SSIM for the slice-based and axial-coupled models.. a — PSNR, b — SSIM
Overall, the choice of \(R\) should be guided by the downstream task. When the analysis relies on weak cellular or nuclear signals, separation of neighboring cells, or accurate 3D/4D tracking, a smaller \(R\) is preferable to preserve fine structural details. When acquisition efficiency is more important and the task can tolerate some loss of fine-scale information, a larger \(R\) may be acceptable. The axial-coupled model does not remove this acquisition-reconstruction trade-off, but it can partially offset the loss of measurement information at stronger compression.
5 Conclusion and future work↩︎
We proposed a PnP-ADMM framework for volumetric reconstruction in CS-LSM. The method handles binary-mask encoded axial measurements and flexibly incorporates any off-the-shelf denoiser, encompassing classical and DL-based methods alike. We developed a slice-based approach as a tractable baseline, while the axial-coupled extension exploits inter-slice continuity to improve reconstruction quality across the volume. Furthermore, we established subsequential convergence of the proposed scheme under a weakly convex regularization assumption, with every accumulation point shown to be a stationary point of the corresponding variational formulation. Experiments on zebrafish-heart data show that the proposed framework can recover cellular structures from compressed measurements while accommodating different denoising methods. The axial-coupled model generally improves reconstruction quality and volumetric continuity over the slice-based model. Future work will investigate tuning-free approaches based on deep unrolling [60] and reinforcement learning [61], as well as extend the proposed framework to 4D spatiotemporal reconstruction. On the theoretical side, establishing convergence guarantees for the proposed axial-coupled PnP framework, beyond the proximal-based denoiser regime, remains an important open problem.
References↩︎
School of Data and Information Sciences, The University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA ().↩︎
Department of Mathematics, The University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA (). Jia and Gong contributed equally to this work.↩︎
Department of Bioengineering, The University of Texas at Dallas, Richardson, TX 75080, USA (, , ).↩︎
Department of Mathematics, School of Data and Information Sciences, The University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA ().↩︎
Submitted to the editors DATE.↩︎
This vectorization is introduced only for notational convenience; in implementation, each slice is kept in its original 2D form.↩︎
The Gaussian noise model is a common approximation that is well-justified at sufficiently high photon counts, which is consistent with our hardware setting.↩︎