Stochastic Multiresolution Image Sketching for Inverse Imaging Problems

TL;DR

Stochastic Multiresolution Image Sketching (ImaSk) addresses the computational bottleneck of high-dimensional inverse imaging by introducing image-domain sketches to form cheaper forward-operator realizations $\mathbf{K}_i=\mathbf{K}\mathbf{S}_i$. It recasts the regularized problem $\min_{\mathbf{x}} \frac{1}{2}\|\mathbf{K}\mathbf{x}-\mathbf{b}\|^2 + R(\mathbf{x})$ as a saddle-point problem and employs a SAGA-like gradient estimator with random multiresolution updates, preserving unbiasedness while reducing per-iteration cost. The authors prove linear convergence for linear forward models under a $\mu$-strongly convex regularizer and demonstrate CT reconstruction gains: increasing the number of resolutions speeds up convergence in wall-clock time. They also introduce ImaSk-Seq, a sequential-update variant with extrapolation that improves early convergence, and show applicability to both smooth ($L_2$) and non-smooth (TV) regularizations. Overall, ImaSk provides a scalable, provably convergent framework for fast, high-resolution image reconstruction with potential extensions to non-linear imaging and measurement-domain sketching.

Abstract

A challenge in high-dimensional inverse problems is developing iterative solvers to find the accurate solution of regularized optimization problems with low computational cost. An important example is computed tomography (CT) where both image and data sizes are large and therefore the forward model is costly to evaluate. Since several years algorithms from stochastic optimization are used for tomographic image reconstruction with great success by subsampling the data. Here we propose a novel way how stochastic optimization can be used to speed up image reconstruction by means of image domain sketching such that at each iteration an image of different resolution is being used. Hence, we coin this algorithm ImaSk. By considering an associated saddle-point problem, we can formulate ImaSk as a gradient-based algorithm where the gradient is approximated in the same spirit as the stochastic average gradient amélioré (SAGA) and uses at each iteration one of these multiresolution operators at random. We prove that ImaSk is linearly converging for linear forward models with strongly convex regularization functions. Numerical simulations on CT show that ImaSk is effective and increasing the number of multiresolution operators reduces the computational time to reach the modeled solution.

Figures (10)

• Figure 1: Illustration of multiresolution sketches. The sketches (here $r=4$) are averaging over small neighborhoods, see \ref{['sec:sketches']} for details. (a) Ground truth image $\mathbf{x}^*$ at full resolution $d = 512^2$, (b, c, d) increasingly lower resolution mappings $\mathbf{S}_i \mathbf{x}^*$ and (e) resolution compensation. All images are shown on the same color range $[-0.1, 1.1]$.
• Figure 1: Visualization of the sinogram decomposition $\mathbf{K}_i\mathbf{x}^*$ with $\mathbf{x}^*$ the ground truth XCAT image in \ref{['fig:CTsketch_ASTRA_Idecim']}(a), $\mathbf{K}$ the parallel X-ray CT operator with $100$ projections and $\sqrt{2d} \approx 724$ detectors implemented with ASTRA Toolbox at image resolutions $\sqrt{d}_i= [64, 128, 256, 512]$, with $i = 1,\ldots,r-1, \, r=4$ and probabilities $p_i=[0.3, 0.3, 0.2], \, p_r = 0.2$: (a) full sinogram $\mathbf{K}\mathbf{x}^*$ from the high resolution image, $\sqrt{d} = 512$, (b-d) sinogram decomposition through low resolution image mappings $\mathbf{K}_i\mathbf{x}^* = \mathbf{R}_i\mathbf{T}_ix^*$ (intensity range $[-0.3, 1.1]$) and (f) error between the full sinogram and the sketched decomposition $\mathbf{K}\mathbf{x}^*-\sum_{i=1}^{r}p_i\mathbf{K}_i\mathbf{x}^*$ (intensity range $[-2, 2]\times 10^{-2}$).
• Figure 1: Quantitative ablation study for ImaSk-Seq for AAPM image. (a) Distance of the ImaSk-Seq estimate $\mathbf{x}^t$ respect to the PDHG solution $\mathbf{x}^{\natural}$, (b) PSNR respect to $\mathbf{x}^*$ in dB.
• Figure 2: Analysis of the computational time of the CT forward operator $\mathbf{R}_i\mathbf{x}$ in ASTRA for different input resolutions $\sqrt{d_i}=2^i,\; i=1,\ldots, 10$ and linear interpolation; mean and standard deviation obtained using 100 Monte--Carlo simulations.
• Figure 2: Comparison of ImaSk, ImaSk-Seq with low-resolution square phantom for $r=8$: (a) distance to optimal solution $\mathbf{x}^{\natural}$ (PDHG), (b) PSNR error respect to $\mathbf{x}^*$ in dB.

Theorems & Definitions (22)

• Theorem 2.1
• Remark 3.1
• Lemma 3.2
• Proof 1
• Remark 3.3
• Theorem 3.4
• Lemma 3.5
• Lemma 3.6
• Proof 2
• Lemma 3.7