Stochastic Reconstruction Technique

• Stochastic reconstruction technique is a computational methodology that employs random multiresolution sketching to solve complex inverse imaging problems.
• It integrates stochastic optimization, variance reduction, and saddle-point reformulation to reduce computational cost and ensure convergence in high-dimensional settings.
• Empirical results in computed tomography validate the method's ability to lower per-iteration cost while achieving linear convergence rates compared to conventional solvers.

Stochastic reconstruction techniques constitute a diverse class of computational methodologies that leverage randomness to solve complex inverse problems, particularly where direct, high-resolution or complete-data solutions are computationally prohibitive or fundamentally underdetermined. Recent advances have integrated stochastic optimization, random sketching, and multiresolution analysis to accelerate regularized iterative solvers in large-scale imaging scenarios, exemplified by computed tomography (CT). The ImaSk algorithm ("Image Sketching"), as presented in (Perelli et al., 13 Dec 2024), strategically combines these principles using random multiresolution image-domain operators to reduce per-iteration cost and maintain convergence guarantees for high-dimensional, regularized reconstruction tasks.

1. Image-Domain Sketching and Multiresolution Operators

ImaSk is predicated on the concept of randomly projecting the original optimization variable (the high-dimensional image) into lower-dimensional subspaces via a set of multiresolution "sketch operators" $\{S_1, S_2, ..., S_r\}$. These operators are designed such that their expectation equals the identity: $\sum_{i=1}^r p_i S_i = I,$ where $p_i$ is the probability of selecting each sketch operator at a given iteration. Each $S_i$ typically arises via block-averaging or other downsampling schemes that map the image to a lower resolution, enabling substantial reductions in the computational complexity of applying the forward operator $K$ (e.g., the Radon transform in CT).

The stochastic process is implemented by selecting, at each solver iteration, a sketch operator $S_{i_k}$ (with probability $p_{i_k}$) and computing using $K_{i_k} = K S_{i_k}$ and its adjoint. This provides an unbiased but noisy estimate of the full forward model, ensuring that iterates are, in expectation, consistent with the original high-resolution problem.

2. Saddle-Point Reformulation for Stochastic Variance-Reduced Primal-Dual Updates

To accommodate the stochastic, nonseparable structure introduced by random sketching, the original regularized least-squares problem,

$\min_x \; \tfrac{1}{2}\|K x - b\|^2 + R(x),$

is reformulated as a convex-concave saddle-point problem: $$\min_{x} \max_{y} \Biggl( \sum_{i=1}^r \langle y, K_i x \rangle - f^*(y) + R(x) \Biggr),$$ where $f(y) = \tfrac{1}{2}\|y - b\|^2$ and $f^*$ is its convex conjugate. Defining $A = (K_1, ..., K_r)$, the problem structure naturally admits the application of stochastic variance-reduced primal-dual updates—specifically, SAGA-type memorization of gradient components—operating on randomly selected multiresolution sketches.

The algorithm maintains a table of "memory variables" for each sketch, updating only the selected entry at each iteration while averaging the contribution from all sketches to preserve unbiasedness. This architecture allows efficient stochastic gradient estimation with desirable variance-reduction properties.

3. Image Sketching Updates and Theoretical Guarantees

The primal and dual variables are updated as follows:

• Let $i_k$ be the randomly chosen sketch at iteration $k$ (according to probabilities $p_i$).
• Update the adjoint memory: $\phi_{i}^{(k+1)} = K_i^{T} y^{(k)}$ if $i = i_k$; else, retain previous value.
• Assemble the stochastic gradient for the primal update:

$$\xi^{(k)} = \phi_{i_k}^{(k+1)} - \phi_{i_k}^{(k)} + \sum_{i=1}^r p_i \phi_i^{(k)}$$

• Update $x$ using the proximal operator of $R$ and a preselected stepsize.
• Similarly, update the dual variable $y$.

The paper affirms that in the case of a linear forward model $K$ and strongly convex regularizer $R$, the algorithm converges linearly toward the global optimum. If $\mu$ is the strong convexity constant and $\sigma$ the stepsize, then for suitably chosen parameters,

$$\mathbb{E}[ \| x^K - x^* \|^2 ] \leq C \cdot \theta^K,$$

with $\theta < 1$ determined by the precise choice of sketch operators and probabilities, and the correlations among the multiresolution operators.

4. Numerical Results in Computed Tomography

The effectiveness of ImaSk is validated numerically on CT image reconstruction tasks. Using a set of downsampling operators $T_i$ to define $S_i = T_i^\top T_i$, the forward operator $K$ is efficiently applied at multiple resolutions (e.g., grid sizes $256 \times 256$ vs. $512 \times 512$). Key empirical observations include:

• Computational cost per iteration decreases commensurately with reduced image resolution.
• Increasing the number of available resolutions ($r$) accelerates convergence in total wall time, and the computational advantage scales according to the complexity of each $K_i$.
• Relative error and PSNR curves plotted vs. "full matrix multiplication" equivalents confirm a substantial time savings as $r$ increases.
• The aggregation property $\sum_i p_i S_i = I$ ensures unbiasing of the reconstruction despite per-iteration information loss due to sketching.

5. Comparison to Other Stochastic Inverse Solvers

Unlike data-domain stochastic or subset methods (e.g., batch SGD, SAGA on measurement indices), ImaSk applies randomness in the image domain. Key comparisons include:

• Lower per-iteration computational cost, as downsampling directly reduces the complexity of the linear projection operations ($K_i x$).
• Built-in variance reduction and linear convergence rates analogous to SAGA, thanks to the memory-augmented stochastic updates.
• Flexibility in trading off accuracy and computation by tuning the set of sketch operators and their selection probabilities.

However, the theoretical guarantees rest on strong convexity assumptions and linear forward models, though empirical evidence with non-strongly convex penalties (like TV) shows favorable performance.

6. Generalizability and Applications

The ImaSk paradigm is not limited to CT but generalizes to any inverse problem with a linear or mildly nonlinear forward model, where forward-map evaluations are expensive:

• PET, MRI, and other large-scale tomographic modalities.
• Inverse problems in remote sensing, industrial NDT, or neuroimaging with large spatial domains.
• Problems requiring regularized optimization with computational constraints, particularly those where a hierarchy of resolutions can be naturally defined.

The stochastic multiresolution update strategy is particularly amenable to parallel and distributed implementations and may be further extended to hybrid approaches (e.g., combining measurement- and image-domain sketching).

7. Future Directions and Potential Extensions

Several extensions of the ImaSk approach are indicated:

• Adapting the saddle-point and variance-reduction structure to nonlinear or nonconvex regularization, including deep learning–based image priors.
• Hybrid schemes that also incorporate stochastic sampling in the data domain.
• Adaptive multiresolution schemes, where the resolution hierarchy and probabilities are adjusted online in response to task-specific criteria.
• Application to problems where model evaluations (e.g., forward PDE solutions) dominate cost, using reduced-order or physics-informed surrogates as sketch operators.

The ImaSk stochastic reconstruction technique, by integrating randomized multiresolution analysis with rigorous optimization theory, contributes a scalable, theoretically principled, and empirically validated solution for large-scale inverse imaging and related high-dimensional reconstruction problems (Perelli et al., 13 Dec 2024).