Volterra example with iterative regularised Gauss-Newton semismooth solver

The discrete Volterra operator. The Volterra operator \(V_n\) is defined as

\[(V_n f)(x) = \int_0^x f(t)^n dt.\]

Its discrete form, using a Riemann sum, is simply

\[(V_n x)_i = h \sum_{j \leq i} x_j^n,\]

where \(h\) is the grid spacing. \(V_1\) is linear.

Therefore, we create a subclass of regpy.Operator type in the next code block. In the following create new operators of that type can simply be done by initiating with that class. This class will have the following parameters. (Side remark you may also do the definition in a separate file for this example check out volterra.py in this example folder. )

---

Parameters

domain : regpy.vecsps.UniformGridFcts > The domain on which the operator is defined. Must be one-dimensional.

exponent : float > The exponent (n). Default is 1.

import numpy as np
import matplotlib.pyplot as plt

from regpy.operators import Operator
from regpy.vecsps import UniformGridFcts


class Volterra(Operator):
    def __init__(self, domain, exponent=1):
        assert isinstance(domain, UniformGridFcts)
        assert domain.ndim == 1
        self.exponent = exponent
        """The exponent."""
        super().__init__(domain, domain, linear=(exponent == 1))

    def _eval(self, x, differentiate=False):
        if differentiate:
            self._factor = self.exponent * x**(self.exponent - 1)
        return self.domain.volume_elem * np.cumsum(x**self.exponent)

    def _derivative(self, x):
        return self.domain.volume_elem * np.cumsum(self._factor * x)

    def _adjoint(self, y):
        x = self.domain.volume_elem * np.flipud(np.cumsum(np.flipud(y)))
        if self.linear:
            return x
        else:
            return self._factor * x

Initialization

Construction of a Volterra operator with exponent \(n=2\) on the interval \([0,2\pi]\) discretized by a UniformGridFcts with 200 points. Note, the operator will be a non-linear operator. Moreover, choose the exact solution

\[f^\dagger(x) = (1-\cos(x))/2\]

for which the exact solution is given as

\[g^\dagger = (3x-4\sin(x)+\cos(x)\sin(x))/8.\]

grid = UniformGridFcts(np.linspace(0, 2 * np.pi, 200))
op = Volterra(grid,exponent=2)

exact_solution = (1-np.cos(grid.coords[0]))**2/4
plt.figure()
plt.plot(grid.coords[0],exact_solution,color="blue",label='exact solution')

exact_data = (3*grid.coords[0] - 4*np.sin(grid.coords[0]) + np.cos(grid.coords[0])*np.sin(grid.coords[0]))/8
noise = 0.3 * op.domain.randn()
data = exact_data + noise
plt.plot(grid.coords[0], exact_data, color="green", label='exact data')
plt.plot(grid.coords[0], data, color="lightgreen", label='noisy data')
plt.legend()

print("Maximal error of operator evaluation on exact solution to exact data: {:.4E}".format(np.max(np.abs(op(exact_solution)-exact_data))))
plt.show()

Maximal error of operator evaluation on exact solution to exact data: 6.3814E-01

Regularization

First we define the regularization setting using the RegularizationSetting and choosing \(L^2\) spaces defined on the UniformGridFcts.

As a next step we define the we choose an initial guess, in this case we choose something that is almost zero.

Use the Semi-Smooth Iterative Regularized Gauss Newton Method implemented in regpy.solvers.nonlinear.irgnm_semismooth to reconstruct the exact solution from the above constructed noisy data data. This method requires an a-priori lower and upper bound on the exact solution here chosen to be \(0\) and \(1.5\). In each step the initial regularization parameter regpar of the solver gets refined by regpar_step. We stop the iterative algorithm after at most \(100\) steps and have as early stopping criterion the discrepancy rule implemented. This can be easily done by summing the two instances of the regpy.stoprules.

After everything is defined run the solver with the specified stopping rule using the method run() of the solver.

from regpy.hilbert import L2
from regpy.solvers import RegularizationSetting
from regpy.solvers.nonlinear.irgnm_semismooth import IrgnmSemiSmooth
import regpy.stoprules as rules

setting = RegularizationSetting(op, L2, L2)

init = op.domain.ones()*0.05

solver = IrgnmSemiSmooth(
    setting, # Regularization setting
    data, # noisy data constructed above
    psi_minus=0, # lower bound
    psi_plus=1.2, # upper bound
    regpar=1, # initial regularization parameter
    regpar_step=0.7, # regularization parameter step
    init=init, # initial guess
)
stoprule = (
    rules.CountIterations(100) +
    rules.Discrepancy(
        setting.h_codomain.norm, data,
        noiselevel=setting.h_codomain.norm(noise),
        tau=1.10
    )
)

recos = []
recos_data = []
for reco, reco_data in solver.while_(stoprule):
    recos.append(1*reco)
    recos_data.append(reco_data)

import matplotlib.animation

plt.rcParams["animation.html"] = "jshtml"
plt.rcParams['figure.dpi'] = 150
plt.ioff()
fig, axs = plt.subplots(1,2)

def animate(t):
    axs[0].cla()
    axs[0].plot(exact_solution,color="blue")
    axs[0].plot(recos[t],color="lightblue")
    axs[0].set_ylim(ymin=0,ymax=1.5)

    axs[1].cla()
    axs[1].plot(exact_data,color="green", label='exact')
    axs[1].plot(recos_data[t],color="orange", label='reco')
    axs[1].plot(data, color="lightgreen", label='measured')
    axs[1].set_ylim(ymin=-0.2,ymax=2.5)

matplotlib.animation.FuncAnimation(fig, animate, frames=len(recos))

2025-08-13 09:04:23,294 INFO CountIterations      :: iteration = 1 / 100
2025-08-13 09:04:23,295 INFO Discrepancy          :: relative discrepancy = 5.14, tolerance = 1.10
2025-08-13 09:04:23,296 INFO IrgnmSemiSmooth      :: Running inner Tikhonov solver.
2025-08-13 09:04:23,299 INFO CountIterations      :: iteration = 2 / 100
2025-08-13 09:04:23,300 INFO Discrepancy          :: relative discrepancy = 1.83, tolerance = 1.10
2025-08-13 09:04:23,300 INFO IrgnmSemiSmooth      :: Running inner Tikhonov solver.
2025-08-13 09:04:23,303 INFO IrgnmSemiSmooth      :: Running inner Tikhonov solver.
2025-08-13 09:04:23,304 INFO IrgnmSemiSmooth      :: Running inner Tikhonov solver.
2025-08-13 09:04:23,307 INFO CountIterations      :: iteration = 3 / 100
2025-08-13 09:04:23,307 INFO Discrepancy          :: relative discrepancy = 1.20, tolerance = 1.10
2025-08-13 09:04:23,308 INFO IrgnmSemiSmooth      :: Running inner Tikhonov solver.
2025-08-13 09:04:23,310 INFO IrgnmSemiSmooth      :: Running inner Tikhonov solver.
2025-08-13 09:04:23,313 INFO CountIterations      :: iteration = 4 / 100
2025-08-13 09:04:23,314 INFO Discrepancy          :: relative discrepancy = 1.02, tolerance = 1.10
2025-08-13 09:04:23,314 INFO CombineRules         :: Rule Discrepancy(noiselevel=0.745453482921017, tau=1.1) triggered.
2025-08-13 09:04:23,315 INFO IrgnmSemiSmooth      :: Solver converged after 3 iteration.
2025-08-13 09:04:23,329 INFO matplotlib.animation :: Animation.save using <class 'matplotlib.animation.HTMLWriter'>

Once Loop Reflect