Deconvolution Problems¶
We consider the two-dimensional deconvolution problems to find a non-negative function f given data
\[d \sim \mathrm{Pois}(h*f)\]
with a non-negative convolution kernel \(h\), and \(\mathrm{Pois}\) denotes the element-wise Poisson distribution.
[1]:
import sys
from pathlib import Path
import numpy as np
from regpy.operators import GaussianBlur
from regpy.solvers import Setting
from regpy.solvers.linear import TikhonovCG
from regpy.functionals import *
from regpy.hilbert import L2
example_dir = "../../../../examples/deconvolution"
Path(example_dir).resolve()
sys.path.insert(0, str(example_dir))
from plots import comparison_plot, convergence_plot, plot_recos
from test_images import mixed
Generate synthetic data with impulsive noise¶
[2]:
grid, exact_sol = mixed(M=256,N=256,c_bubbles=1.,fac=50.)
r"""grid is the underlying UniformGridFcts vector space, and exact_sol the exact solution."""
a=0.05
conv = GaussianBlur(grid,a,pad_amount=16)
r"""Convolution operator \(f\mapsto h*f\) for the convolution kernel \(h(x)=\exp(-|x|_2^2/a^2)\)."""
blur = conv(exact_sol)
blur[blur<0] = 0.
data = blur.copy()
n,m=grid.shape
for j in range(64**2):
nn = np.random.randint(n)
mm = np.random.randint(m)
data[nn,mm] = np.random.randint(2000)-1000
comparison_plot(grid,exact_sol,data,title_left='noisy measurement data')
Reconstruction with a quadratic data fidelity term¶
… fails completely for this data!
[3]:
alpha = 1e-4
setting = Setting(conv,L2,L2)
solver = TikhonovCG(setting=setting, data = data, regpar = alpha,reltolx=1e-6)
fal,_ = solver.run()
comparison_plot(grid,exact_sol,fal,title_left='quadratic Tikhonov regularization')
2025-12-19 19:09:07,703 WARNING Setting :: Setting does not contain any explicit data.
2025-12-19 19:09:10,932 INFO TikhonovCG :: it.745 kappa=25.918371183841256 err/Tol rel X:9.8e-07/1.0e-06
2025-12-19 19:09:10,933 INFO TikhonovCG :: Solver converged after 745 iteration.
setting with Huber data fidelity and quadratic penalty term¶
We now set up a generalized Tikhonov regularization setting
\[\hat{f} \in \mathrm{argmin}_f\left[\frac{1}{\alpha}H_{\sigma}(Tf-g^{\mathrm{obs}})+\|f\|^2\right]\]
where the Huber data fidelity term \(H_{\sigma}\) is quadratic on coordinates in \([-\sigma,\sigma]\) and \(|\cdot|-\sigma/2\) for large coordinates.
[4]:
from regpy.functionals import Huber,LppPower
sigma = 0.1
Sdat = (1./sigma)*Huber(grid,sigma=sigma,eps =1e-10)
L2pen = LppPower(grid,p=2)
huber_setting = Setting(conv,L2pen,Sdat,alpha,data=data)
display the available methods
[5]:
huber_setting.display_all_methods(full_names=False)
Applicable methods:
FB ForwardBackwardSplitting used with step length tau=1.011e-01.
Expected linear convergence rate: 1.000e+00. linear rate: 1.000e+00
dual_FB ForwardBackwardSplitting used with step length tau=1.011e+00.
Expected linear convergence rate: 1.000e+00. linear rate: 1.000e+00
FISTA FISTA used with convexity parameters mu_R=1.000e-04, mu_S=0.000e+00 and step length tau=1.011e-01.
Expected linear convergence rate: 9.968e-01.
linear rate: 9.968e-01
dual_FISTA FISTA used with convexity parameters mu_R=1.000e-05, mu_S=0.000e+00 and step length tau=1.011e+00.
Expected linear convergence rate: 9.968e-01.
linear rate: 9.968e-01
PDHG Using accelerated version 2 with convexity parameters mu_R=1.000e+00, mu_S*=1.000e-05 and ||T||=9.945e-01.
Expected linear convergence rate: 9.937e-01 linear rate: 9.937e-01
dual_PDHG Using accelerated version 2 with convexity parameters mu_R=1.000e-09, mu_S*=1.000e+04 and ||T||=9.945e-01.
Expected linear convergence rate: 9.937e-01 linear rate: 9.937e-01
ADMM Ergodic rate O(1/n). linear rate: -1.000e+00
dual_SSNewton linear rate: nan
Non-applicable methods:
SSNewton Data functional not quadratic.
[6]:
methods = ['FISTA','dual_FISTA','PDHG','dual_PDHG','dual_FB']
recos = {}
for method in methods:
#huber_setting.set_stopping_rule(method, OptimalityCondStopping(huber_setting,tol=0.01,max_iter=500,logging_level=logging.INFO))
recos[method]= huber_setting.run(method)[0]
2025-12-19 19:09:11,862 INFO FISTA :: FISTA used with convexity parameters mu_R=1.000e-04, mu_S=0.000e+00 and step length tau=1.005e-01.
Expected linear convergence rate: 9.968e-01.
2025-12-19 19:09:17,363 INFO FISTA :: Solver converged after 1000 iteration.
2025-12-19 19:09:17,554 INFO FISTA :: FISTA used with convexity parameters mu_R=1.000e-05, mu_S=0.000e+00 and step length tau=1.005e+00.
Expected linear convergence rate: 9.968e-01.
2025-12-19 19:09:20,247 INFO FISTA :: Solver converged after 485 iteration.
2025-12-19 19:09:20,428 INFO PDHG :: Using accelerated version 2 with convexity parameters mu_R=1.000e+00, mu_S*=1.000e-05 and ||T||=9.936e-01.
Expected linear convergence rate: 9.937e-01
2025-12-19 19:09:27,333 INFO PDHG :: Solver converged after 1000 iteration.
2025-12-19 19:09:27,476 INFO PDHG :: Using accelerated version 2 with convexity parameters mu_R=1.000e-09, mu_S*=1.000e+04 and ||T||=9.948e-01.
Expected linear convergence rate: 9.937e-01
2025-12-19 19:09:29,512 INFO PDHG :: Solver converged after 303 iteration.
2025-12-19 19:09:35,072 INFO ForwardBackwardSplitting :: Solver converged after 1000 iteration.
plot the convergence histories
[7]:
convergence_plot(huber_setting,methods,'Huber data fid., quadratic penalty')
visualize the reconstructions
[8]:
plot_recos(huber_setting,methods,recos,truth=exact_sol)
