Inverse potential problem¶
For the inverse potential the goal is to find a homogeneous source given measurements at a circle enclosing the source. The Operator that maps the shape of a homogeneous heat source to the heat flux measured at some circle outside of the object is defined by the following partial differential equation. The heat distributions satisfies
where \(\partial\Omega\) is the measurement circle and \(K\) is the heat source. The forward operator maps the shape of the heat source to the Neumann data:
References:
Hettlich & W. Rundell “Iterative methods for the reconstruction of an inverse potential problem”, Inverse Problems, 12 (1996) 251–266.
Hohage “Logarithmic convergence rates of the iteratively regularized Gauss–Newton method for an inverse potential and an inverse scattering problem”, Inverse Problems, 13 (1997) 1279–1299.
Import auxiliary libraries and define logging¶
We use numpy arrays as the main data structure for the domain and rely on matplotlib.pyplot as the tool for plots.
To observe the information of the solver while computing we define the logging string and set the log level to info.
[1]:
import logging
import matplotlib.pyplot as plt
import numpy as np
logging.basicConfig(
level=logging.INFO,
format='%(asctime)s %(levelname)s %(name)-40s :: %(message)s'
)
Defining the Operator class¶
We define a class for the potential operator that operators on a uniform grid to discretize the measurements on the circle and relies on the StarTrigDiscr class from regpy as a discretization of star shaped sources.
[2]:
from regpy.operators import Operator
from regpy.vecsps.curve import StarTrigDiscr
from regpy.vecsps import UniformGridFcts
class Potential(Operator):
r"""Operator that maps the shape of a homogeneous heat source to the heat flux measured at some
circle outside of the object. The heat distributions satisfies
.. math::
\begin{cases}
\Delta u = 1_K & \text{ in } \Omega \\
u = 0 & \text{ on } \partial\Omega
\end{cases}
where \(\partial\Omega\) is the measurement circle and \(K\) is the heat source. The operator
maps the shape of the heat source to the Neumann data:
.. math::
\partial K \mapsto \left.\frac{\partial u}{\partial\nu}\right|_{\partial\Omega}.
Attributes
----------
radius : float
The radius of the measurement circle.
N_means : int
The number of equispaced measurement points on the circle.
nforward : int, optional
The order of the Fourier expansion in the forward solver.
Raises
------
ValueError
Will be raised on evaluating the operator if the object radius is negative or penetrates
the measurement circle.
References
----------
- F. Hettlich & W. Rundell "Iterative methods for the reconstruction of an inverse potential
problem", Inverse Problems, 12 (1996) 251–266.
- T. Hohage "Logarithmic convergence rates of the iteratively regularized
Gauss–Newton method for an inverse potential and an inverse scattering problem", Inverse
Problems, 13 (1997) 1279–1299.
"""
def __init__(self, radius, nforward=128, N_means=128, N_ieq=128):
self.radius = radius
"""The measurement radius."""
self.nforward = nforward
"""The Fourier order of the forward solver."""
super().__init__(
domain=StarTrigDiscr(2*N_ieq),
codomain=UniformGridFcts(np.linspace(0, 2*np.pi, N_means, endpoint=False), dtype=complex)
)
k = 1 + np.arange(self.nforward)
k_t = np.outer(k, np.linspace(0, 2*np.pi, self.nforward, endpoint=False))
k_tfl = np.outer(k, self.codomain.coords[0])
self.cosin = np.cos(k_t)
self.sinus = np.sin(k_t)
self.cos_fl = np.cos(k_tfl)
self.sin_fl = np.sin(k_tfl)
def _eval(self, x, differentiate=False, adjoint_derivative=False):
nfwd = self.nforward
self._bd = self.domain.eval_curve(x, nvals=nfwd)
q = self._bd.radius[0]
if q.max() >= self.radius:
raise ValueError('Object penetrates measurement circle')
if q.min() <= 0:
raise ValueError('Radial function negative')
qq = q**2
flux = 1 / (2 * self.radius * nfwd) * np.sum(qq) * self.codomain.ones()
fac = 2 / (nfwd * self.radius)
for j in range(0, (nfwd - 1) // 2):
fac /= self.radius
qq *= q
flux += (
(fac / (j + 3)) * self.cos_fl[j, :] * np.sum(qq * self.cosin[j, :]) +
(fac / (j + 3)) * self.sin_fl[j, :] * np.sum(qq * self.sinus[j, :])
)
if nfwd % 2 == 0:
fac /= self.radius
qq *= q
flux += fac * self.cos_fl[:, nfwd // 2] * np.sum(qq * self.cosin[nfwd // 2, :])
return flux
def _derivative(self, h):
nfwd = self.nforward
q = self._bd.radius[0]
qqh = q * self._bd.derivative(h)
der = 1 / (self.radius * nfwd) * np.sum(qqh) * self.codomain.ones()
fac = 2 / (nfwd * self.radius)
for j in range((nfwd - 1) // 2):
fac /= self.radius
qqh *= q
der += fac * (
self.cos_fl[j, :] * np.sum(qqh * self.cosin[j, :]) +
self.sin_fl[j, :] * np.sum(qqh * self.sinus[j, :])
)
if nfwd % 2 == 0:
fac /= self.radius
qqh *= q
der += fac * self.cos_fl[nfwd // 2, :] * np.sum(qqh * self.cosin[nfwd // 2, :])
return der
def _adjoint(self, g):
nfwd = self.nforward
q = self._bd.radius[0]
qq = q.copy()
adj = 1 / (self.radius * nfwd) * np.sum(g) * qq
fac = 2 / (nfwd * self.radius)
for j in range((nfwd - 1) // 2):
fac /= self.radius
qq *= q
adj += fac * (
np.sum(g * self.cos_fl[j, :]) * (self.cosin[j, :] * qq) +
np.sum(g * self.sin_fl[j, :]) * (self.sinus[j, :] * qq)
)
if nfwd % 2 == 0:
fac /= self.radius
qq *= q
adj += fac * np.sum(g * self.cos_fl[nfwd // 2, :]) * (self.cosin[nfwd // 2, :] * qq)
return self._bd.adjoint(adj.real)
Define Operator and regularization setting¶
Using the above class create the forward operator.
The regularization setting is defined using a Sobolev penalty and \(L^2\) data-fidelity. Note, that the construction here uses the abstract Sobolev which by default is the Sobolev space of with smoothness index 1.
[3]:
from regpy.solvers import RegularizationSetting
from regpy.hilbert import L2, Sobolev
#Forward operator
op = Potential(
radius=1.3
)
setting = RegularizationSetting(op=op, penalty=Sobolev, data_fid=L2)
The exact solution and data¶
For the exact solution we define an star shaped domain with boundary given by $ \sqrt{3\cos{t}^2+1}/2$. We then define the exact data and add noise with 3% relative noise level.
[4]:
#Exact data and Poission data
exact_solution = op.domain.sample(lambda t: np.sqrt(6*np.cos(1.5*t)**2+1)/3)
exact_data = op(exact_solution)
noise = op.codomain.randn()
noise = 0.03*setting.h_codomain.norm(exact_data)/setting.h_codomain.norm(noise) * noise
data = exact_data + noise
fig = plt.figure(figsize=(7,7))
plt.plot(*op.domain.eval_curve(exact_solution).curve[0])
ax = plt.gca()
ax.set_xlim([-1.1, 1.1])
ax.set_ylim([-1.1, 1.1])
plt.show()
Solving the problem¶
With an initial guess that is the circle of radius 1 we use a Newton CG solver to get the solution. For a reference we plot the initial guess with the true objective.
As stopping rules we rely on the discrepancy principle and to prevent an extensively long run we combine it with an iteration maximal count of 100.
[5]:
#Initial guess
init = op.domain.sample(lambda t: 0.2)
fig = plt.figure(figsize=(7,7))
plt.plot(*op.domain.eval_curve(exact_solution).curve[0])
plt.plot(*op.domain.eval_curve(init).curve[0])
ax = plt.gca()
ax.set_xlim([-1.1, 1.1])
ax.set_ylim([-1.1, 1.1])
plt.show()
from regpy.stoprules import CountIterations, Discrepancy
from regpy.solvers.nonlinear.newton import NewtonCG
solver = NewtonCG(
setting,
data,
init = init,
cgmaxit=50,
rho=0.6
)
stoprule = (
CountIterations(100) +
Discrepancy(
setting.h_codomain.norm, data,
noiselevel = setting.h_codomain.norm(noise),
tau=1.1
)
)
[6]:
#Plot function
for reco,reco_data in solver.until(stoprule=stoprule):
fig, axs = plt.subplots(1, 2,figsize=(15,7))
axs[0].set_title('Obstacle')
axs[0].set_xlim([-1.1, 1.1])
axs[0].set_ylim([-1.1, 1.1])
axs[1].set_title('Heat flux')
axs[0].plot(*op.domain.eval_curve(exact_solution).curve[0])
axs[0].plot(*op.domain.eval_curve(reco).curve[0])
axs[1].plot(exact_data.real, label='exact')
axs[1].plot(reco_data.real, label='reco')
axs[1].plot(data.real, label='measured')
axs[1].legend()
axs[1].set_ylim(ymin=0)
plt.show()
2025-08-13 09:03:48,329 INFO NewtonCG :: Inner CG iteration required 1 steps.
2025-08-13 09:03:48,472 INFO CountIterations :: iteration = 1 / 100
2025-08-13 09:03:48,475 INFO Discrepancy :: relative discrepancy = 76.29, tolerance = 1.10
2025-08-13 09:03:48,480 INFO NewtonCG :: Inner CG iteration required 1 steps.
2025-08-13 09:03:48,615 INFO CountIterations :: iteration = 2 / 100
2025-08-13 09:03:48,616 INFO Discrepancy :: relative discrepancy = 13.51, tolerance = 1.10
2025-08-13 09:03:48,620 INFO NewtonCG :: Inner CG iteration required 1 steps.
2025-08-13 09:03:48,812 INFO CountIterations :: iteration = 3 / 100
2025-08-13 09:03:48,813 INFO Discrepancy :: relative discrepancy = 2.75, tolerance = 1.10
2025-08-13 09:03:48,829 INFO NewtonCG :: Inner CG iteration required 8 steps.
2025-08-13 09:03:48,970 INFO CountIterations :: iteration = 4 / 100
2025-08-13 09:03:48,971 INFO Discrepancy :: relative discrepancy = 2.65, tolerance = 1.10
2025-08-13 09:03:48,974 INFO NewtonCG :: Inner CG iteration required 1 steps.
2025-08-13 09:03:49,117 INFO CountIterations :: iteration = 5 / 100
2025-08-13 09:03:49,118 INFO Discrepancy :: relative discrepancy = 1.03, tolerance = 1.10
2025-08-13 09:03:49,118 INFO CombineRules :: Rule Discrepancy(noiselevel=0.012891189621064086, tau=1.1) triggered.
2025-08-13 09:03:49,119 INFO NewtonCG :: Solver converged after 5 iteration.