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')
2026-08-24 15:42:27,404 WARNING Setting :: Setting does not contain any explicit data.
WARNING:Setting:Setting does not contain any explicit data.
2026-08-24 15:42:27,406 WARNING Setting :: Overwriting existing data in setting!
WARNING:Setting:Overwriting existing data in setting!
2026-08-24 15:42:31,132 INFO CombineRules :: it. 746>=1000 | (rel X:9.8e-07>=1.0e-06 & kappa:2.1e+01[All rules triggered.])
Rule AndCombineRules([RelTolXStop(1e-06), KappaTrackAndMinIt(0)]) triggered.
INFO:CombineRules:it. 746>=1000 | (rel X:9.8e-07>=1.0e-06 & kappa:2.1e+01[All rules triggered.])
Rule AndCombineRules([RelTolXStop(1e-06), KappaTrackAndMinIt(0)]) triggered.
2026-08-24 15:42:31,137 INFO TikhonovCG :: Solver converged after 746 iteration.
INFO:TikhonovCG:Solver converged after 746 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.014e-01.
Expected linear convergence rate: 1.000e+00. linear rate: 1.000e+00
dual_FB ForwardBackwardSplitting used with step length tau=1.014e+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.014e-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.014e+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.933e-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.933e-01.
Expected linear convergence rate: 9.937e-01 linear rate: 9.937e-01
ADMM Ergodic rate O(1/n). linear rate: -1.000e+00
SSNewton linear rate: nan
dual_SSNewton linear rate: nan
Non-applicable methods:
TikhCG Penalty term is not convex quadratic.Data fidelity term is not convex quadratic.
dual_TikhCG Penalty term is not convex quadratic.Data fidelity term is not convex 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]
2026-08-24 15:42:31,994 INFO FISTA :: FISTA used with convexity parameters mu_R=1.000e-04, mu_S=0.000e+00 and step length tau=1.010e-01.
Expected linear convergence rate: 9.968e-01.
INFO:FISTA:FISTA used with convexity parameters mu_R=1.000e-04, mu_S=0.000e+00 and step length tau=1.010e-01.
Expected linear convergence rate: 9.968e-01.
2026-08-24 15:42:39,211 INFO CombineRules :: duality gap:6.53e+05 <= 1.0e-01 | it. 1000>=1000
Rule CountIterations(max_iterations=1000) triggered.
INFO:CombineRules:duality gap:6.53e+05 <= 1.0e-01 | it. 1000>=1000
Rule CountIterations(max_iterations=1000) triggered.
2026-08-24 15:42:39,212 INFO FISTA :: Solver converged after 1000 iteration.
INFO:FISTA:Solver converged after 1000 iteration.
2026-08-24 15:42:39,374 INFO FISTA :: FISTA used with convexity parameters mu_R=1.000e-05, mu_S=0.000e+00 and step length tau=1.015e+00.
Expected linear convergence rate: 9.968e-01.
INFO:FISTA:FISTA used with convexity parameters mu_R=1.000e-05, mu_S=0.000e+00 and step length tau=1.015e+00.
Expected linear convergence rate: 9.968e-01.
2026-08-24 15:42:42,682 INFO CombineRules :: duality gap:9.87e-02 <= 1.0e-01 | it. 478>=1000
INFO:CombineRules:duality gap:9.87e-02 <= 1.0e-01 | it. 478>=1000
2026-08-24 15:42:42,684 INFO FISTA :: Solver converged after 478 iteration.
INFO:FISTA:Solver converged after 478 iteration.
2026-08-24 15:42:42,846 INFO PDHG :: Using accelerated version 2 with convexity parameters mu_R=1.000e+00, mu_S*=1.000e-05 and ||T||=9.965e-01.
Expected linear convergence rate: 9.937e-01
INFO:PDHG:Using accelerated version 2 with convexity parameters mu_R=1.000e+00, mu_S*=1.000e-05 and ||T||=9.965e-01.
Expected linear convergence rate: 9.937e-01
2026-08-24 15:42:51,935 INFO CombineRules :: duality gap:2.86e+00 <= 1.0e-01 | it. 1000>=1000
Rule CountIterations(max_iterations=1000) triggered.
INFO:CombineRules:duality gap:2.86e+00 <= 1.0e-01 | it. 1000>=1000
Rule CountIterations(max_iterations=1000) triggered.
2026-08-24 15:42:51,937 INFO PDHG :: Solver converged after 1000 iteration.
INFO:PDHG:Solver converged after 1000 iteration.
2026-08-24 15:42:52,140 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
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
2026-08-24 15:42:55,095 INFO CombineRules :: duality gap:9.95e-02 <= 1.0e-01 | it. 303>=1000
INFO:CombineRules:duality gap:9.95e-02 <= 1.0e-01 | it. 303>=1000
2026-08-24 15:42:55,097 INFO PDHG :: Solver converged after 303 iteration.
INFO:PDHG:Solver converged after 303 iteration.
2026-08-24 15:43:03,602 INFO CombineRules :: duality gap:2.43e-01 <= 1.0e-01 | it. 1000>=1000
Rule CountIterations(max_iterations=1000) triggered.
INFO:CombineRules:duality gap:2.43e-01 <= 1.0e-01 | it. 1000>=1000
Rule CountIterations(max_iterations=1000) triggered.
2026-08-24 15:43:03,604 INFO ForwardBackwardSplitting :: Solver converged after 1000 iteration.
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)
