Learned Regularization for Inverse Problems: Insights from a Spectral Model

In this respository, we document the numerical examples from our paper Learned Regularization for Inverse Problems: Insights from a Spectral Model.

Content

Data

Requirements

Experiments

References

Citation

Data

In order to run the code smoothly and without too many package dependencies, some pre-saved data is required. The respective folder spectralData can be downloaded here:
https://syncandshare.desy.de/index.php/s/rSysTmdGAtYAf7r

The folder has a size of approximately 5.6 GB and contains:

• Three operators that implement the Radon transform for inputs of different dimensions, created with the radon package, see also [1]. To simplify the computations we perform in our experiments, each operator contains the matrices $U,V$ and the vector $\sigma$ that form the singular value decomposition
  $$A = V\mathrm{diag}(\sigma)U^{T}.$$
  The pre-saved operators have the following specifications:
  • radon_32.obj (21.6 MB) with forward operator $A: \mathbb{R}^{32\times 32} \rightarrow \mathbb{R}^{32\times 47}$
  • radon_64.obj (343 MB) with forward operator $A: \mathbb{R}^{64\times 64} \rightarrow \mathbb{R}^{64\times 93}$
  • radon_128.obj (5.3 GB) with forward operator $A: \mathbb{R}^{128\times 128} \rightarrow \mathbb{R}^{128\times 183}$

from radon import radon_matrix 
from ops import svd_op
import torch

res = 16                     # your favorite resolution

A = radon_matrix(img = torch.zeros((res,res)), 
                thetas = torch.linspace(0.0, torch.pi, res)[:-1])

op16 = svd_op(A,res)     # computes the SVD and implements forward, inverse and adjoint for 2-D inputs

• Three lists (pi32.txt, pi64.txt, pi128.txt) of coefficients $\{\Pi_n^*\}_{n \in \mathbb{N}}$ corresponding to the three pre-saved operators. The data distribution $\pi^*$ is described by the first 2400 images of the LoDoPaB CT data set, see also [2].
  The code we have used to compute the coefficients is documented in coeffs.py and can be applied to any other combination of data and operator.
• An example image test_img.png, taken from the LoDoPaB CT data set.

Requirements

Using the pre-saved data, running the code only requires the packages matplotlib, numpy, and scikit-image.
You can install these packages by running

pip install -r requirements.txt

in the project directory.

Experiments

The aim of the experiments we provide here, is to compare different data-driven spectral reconstruction operators to solve inverse problems. More precisely, for an inverse problem $y= Ax + \epsilon$, we consider reconstruction operators

$$ R_{\mu}(y) = \sum_{n = 1}^N g_n(\mu) \langle y, v_n \rangle u_n = U\mathrm{diag}(g(\mu))V^T y.$$

Here, $A = V\mathrm{diag}(\sigma)U^{T}$ is the singular value decompositon of $A$. The n-th columns of $U$ and $V$ are denoted by $u_n$ and $v_n$.

The filter $g(\mu) = \{g_n(\mu)\}_{n \in \mathbb{N}}$ is now optimized with respect to three different learning paradigms (which we call mse, adv and post in the code),
with training noise $\epsilon$ drawn from the distribution $\mu$ and data $x$ drawn from the distribution $\pi^*$.
For all considered approaches, the optimal filter coefficients are of the form

$$ g_n(\mu) = \frac{\Pi^*_n\sigma_n}{\Pi_n^*\sigma_n^2 + \lambda_n(\sigma_n, \Pi^*_n, \Delta_n(\mu))}.$$

Here, $\sigma_n$ denotes again the singular values, and the data-driven, or, respectively, noise-driven coefficients are given by

$$ \Pi_n^* = \mathbb{E}_{x \sim \pi^*}[\langle x, u_n\rangle^2] \quad \text{and}\quad \Delta_n(\mu) = \mathbb{E}_{\epsilon \sim \mu}[\langle \epsilon, v_n\rangle^2].$$

The specific coefficients for the different approaches are then given by

$$ \color{darkorange}{\lambda_n^{\text{mse}}(\mu)} = \Delta_n(\mu),\qquad \color{teal}{\lambda_n^{\text{adv},3/8}(\mu)}= \frac{\Delta_n(\mu)}{3\sigma_n^2\Pi_n^* + \Delta_n(\mu)}, \qquad \color{red}{\lambda_n^{\text{post}(\mu)}} = \sigma_n^2 \Delta_n(\mu).$$

Our driving questions are:

1. Are the data-driven reconstruction operators actually regularizers, i.e., stable?
2. If yes: Is the regularization convergent, i.e., do we converge to a "ground truth"-reconstruction as the noise level tends to zero?
3. If yes: Is the regularization still convergent if we test it on noise drawn from $\nu \neq \mu$?

In the experiments, we perform CT-reconstruction, where the forward operator $A$ is given by the Radon transform.

The main parts of this repository are the two Jupyter Notebooks continuity.ipynb and convergence.ipynb.

• In continuity.ipynb (which complements Section 6.1 of our paper) we document the experiments on the reconstruction error for a fixed and known noise model on three different image resolutions.
• In convergence.ipynb (which complements Section 6.2 of our paper) we document the experiments on the reconstruction error for a fixed image resolution but three different noise models and decaying noise level.

References

[1]
Kabri, S., Auras, A., Riccio, D., Bauermeister, H., Benning, M., Moeller, M., Burger, M.
Convergent Data-Driven Regularizations for CT Reconstruction. Commun. Appl. Math. Comput. (2024).
[2]
Leuschner, J., Schmidt, M., Otero Baguer, D., Maass, P.
LoDoPaB-CT, a benchmark dataset for low-dose computed tomography reconstruction.
Scientific Data, 8(109). (2021).

Citation

If you would like to cite our book chapter, you can copy the following bib entry:

@inbook{Burger2025learned,
url = {https://doi.org/10.1515/9783111251233-002},
title = {Learned regularization for inverse problems},
booktitle = {Data-driven Models in Inverse Problems},
author = {Martin Burger and Samira Kabri},
editor = {Tatiana A. Bubba},
publisher = {De Gruyter},
address = {Berlin, Boston},
pages = {39--72},
doi = {doi:10.1515/9783111251233-002},
isbn = {9783111251233},
year = {2025},
lastchecked = {2025-05-20}
}