Operators in RegPy

This tutorial expands on the usage guide and explains how to define custom operators using regpy.operators.Operator.

In RegPy, an operator represents a possibly nonlinear mapping between vector spaces \(\mathbb{X}\) and \(\mathbb{Y}\):

\[F\colon \mathbb{X}\to\mathbb{Y}.\]

The spaces \(\mathbb{X}\) and \(\mathbb{Y}\) are instances of regpy.vecsps.VectorSpaceBase. Their vectors are duck typed and need not be NumPy arrays; see Spaces. Complex coefficient spaces are treated as real vector spaces, so a finite-dimensional domain and codomain can be identified with \(\mathbb{R}^N\) and \(\mathbb{R}^M\), respectively.

A linear operator \(T\colon \mathbb{X}\to\mathbb{Y}\) can be represented by a matrix \(\underline{T}\in\mathbb{R}^{M\times N}\). Forming this matrix may be inefficient or impossible, whereas applying it to a vector is often feasible. RegPy represents this action through _eval and the action of the standard vector-space adjoint through _adjoint.

Equivalently, _adjoint represents the dual operator \(T'\colon\mathbb{Y}'\to\mathbb{X}'\) after identifying each space with its dual through the standard real pairing:

\[\langle u,v\rangle = \operatorname{Re}(\operatorname{vdot}(u,v)).\]

The same real-pairing convention is used throughout RegPy, including for functionals.

Thus _eval and _adjoint must satisfy

T.codomain.vdot(T(x), y).real == T.domain.vdot(x, T.adjoint(y)).real

If \(T\) is \(\mathbb{C}\)-linear and represented by a matrix \(\underline{T}\in\mathbb{C}^{M\times N}\), then T.adjoint is represented by the conjugate transpose of \(\underline{T}\). In that case, the corresponding complex identity also holds before taking real parts.

Regularization methods often equip the spaces with additional Hilbert-space structures. If their inner products are represented by the Gram matrices \(G_{\mathbb{X}}\) and \(G_{\mathbb{Y}}\), then the corresponding Hilbert-space adjoint is

\[T^{\ast} = G_{\mathbb{X}}^{-1} T' G_{\mathbb{Y}},\]

where \(T'\) denotes the adjoint with respect to the standard vector-space pairings.

This decomposition motivates RegPy’s design: an operator adjoint is defined with respect to the standard pairings supplied by its vector spaces. Additional Hilbert-space geometry is introduced separately through Gram operators. The operator representation therefore remains independent of the inner products, data-fidelity terms, and penalties used by a particular inverse problem.

Using existing operators

The easiest way to construct an operator is often to combine existing classes from regpy.operators. The module includes multiplication, Fourier transform, convolution, direct-sum, and composition operators.

op_1 = Some_Operator(...)
op_2 = Some_other_Operator(...)

my_op = op_1 * op_2  # requires op_2.codomain == op_1.domain
shifted_op = my_op + my_op.codomain.rand()

Common operator operations include:

• a * op1 + b * op2: linear combination;

• op1 * op2: composition;

• op * factor: composition with pointwise multiplication in the domain;

• op + offset: shift in the codomain.

Linear operators

Every operator requires domain and codomain instances derived from regpy.vecsps.VectorSpaceBase. They may be supplied to the constructor or constructed from its parameters. These objects describe vector membership and storage operations; they do not prescribe a Hilbert-space inner product.

The initialization

For example, an operator between two uniform two-dimensional grids may accept axis specifications (start, end, number) and construct regpy.vecsps.UniformGridFcts instances:

def __init__(self, d_1, d_2, cd_1, cd_2):
    domain = UniformGridFcts(d_1, d_2)
    codomain = UniformGridFcts(cd_1, cd_2)
    super().__init__(domain=domain, codomain=codomain, linear=True)

For a linear operator, you need to implement two methods:

• _eval computes the action of the forward operator.

• _adjoint computes the action of its adjoint.

Example

Consider a two-dimensional Fourier transform on a centered square grid. Let d = (-1, 1, 100) specify each domain axis. The dual grid spacing determines the codomain:

def __init__(self, d):
    domain = UniformGridFcts(d, d, dtype=complex)
    cd = (-1 / (2 * domain.spacing[0]),
          1 / (2 * domain.spacing[0]),
          domain.shape[0])
    codomain = UniformGridFcts(cd, cd, dtype=complex)
    super().__init__(
        domain=domain,
        codomain=codomain,
        linear=True,
    )

The evaluation method

For a linear operator, _eval receives a single vector x. Users invoke the public interface as op(x); regpy.operators.Operator checks that x belongs to the domain and that the returned vector belongs to the codomain. An implementation of _eval may therefore assume a valid input, but it must return a freely modifiable codomain vector.

Example

For the two-dimensional Fourier transform, the NumPy implementation is:

def _eval(self, x):
    return np.fft.fftshift(np.fft.fftn(np.fft.ifftshift(x), norm='ortho'))

The adjoint evaluation method

The _adjoint method follows the same convention. Users call op.adjoint(y); the public interface validates y against the codomain and validates the result against the domain. The adjoint is always taken with respect to the standard real pairings of the two vector spaces.

As mentioned above, RegPy defines adjoints with respect to the standard real pairing. For vectors with complex coefficients, this pairing is

\[\langle x,y\rangle_{\mathbb{R}} = \operatorname{Re}(x^*y).\]

If _eval is represented by a matrix \(\underline{T}\), the following cases are useful:

• If the domain and codomain have real coefficients, _adjoint implements \(y\mapsto \underline{T}^{\top}y\).

• If the domain and codomain have complex coefficients and _eval is complex-linear, _adjoint implements \(y\mapsto \underline{T}^{*}y\).

• If the domain has real coefficients, the codomain has complex coefficients, and _eval implements \(x\mapsto \underline{T}x\), then _adjoint implements \(y\mapsto \operatorname{Re}(\underline{T}^{*}y)\).

• If the domain has complex coefficients, the codomain has real coefficients, and _eval implements \(x\mapsto \operatorname{Re}(\underline{T}x)\), then _adjoint implements \(y\mapsto \underline{T}^{*}y\).

We recommend checking every new _adjoint implementation with regpy.util.operator_tests.test_adjoint().

For complex coefficient spaces, let \(\underline{T}\) be the complex matrix of a complex-linear operator and let \(G_{\mathbb{X}}\) and \(G_{\mathbb{Y}}\) be Hermitian positive-definite Gram matrices. The Hilbert-space adjoint is then represented by

\[T_{\mathbb{X},\mathbb{Y}}^{\ast} = G_{\mathbb{X}}^{-1}\underline{T}^{*}G_{\mathbb{Y}}.\]

This is consistent with the usual complex Hilbert-space adjoint and with RegPy’s real-pairing convention.

Example

For the example above, the inverse Fourier transform defines the adjoint:

def _adjoint(self, y):
    return np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(y), norm='ortho'))

Allocating and in-place implementations

Each fundamental operator action has an allocating and an in-place variant:

• _eval(x, ...) and _ieval(x, out, ...) evaluate the operator;

• _derivative(h, ...) and _iderivative(h, out, ...) apply the derivative at the current linearization point;

• _adjoint(y, ...) and _iadjoint(y, out, ...) apply the corresponding adjoint.

The leading i means that the result is written into the supplied out vector. These are private implementation methods; callers continue to use the public operator interface and pass an output buffer with the out keyword:

y = op.codomain.empty()
op(x, out=y)

x_adjoint = op.domain.empty()
op.adjoint(y, out=x_adjoint)

For a nonlinear operator, the derivative returned by linearize has the same public interface:

value, derivative = op.linearize(x)

derivative_value = op.codomain.empty()
derivative(h, out=derivative_value)

adjoint_value = op.domain.empty()
derivative.adjoint(y, out=adjoint_value)

An operator normally implements only one method from each allocating/in-place pair. The base class derives the other variant automatically. If only the in-place method is implemented, the allocating variant creates a vector in the appropriate output space and passes it as out. If only the allocating method is implemented, the default in-place variant evaluates it and transfers the result into out using the output vector space’s copyto method.

Implementing an in-place variant is useful when an algorithm can write directly to existing storage and thereby avoid an allocation or an additional copy. Its contract is:

• out belongs to the method’s output space: the codomain for _ieval and _iderivative, and the domain for _iadjoint;

• the method completely overwrites out and does not rely on its previous contents;

• it returns the exact same object supplied as out.

Whenever the domain and codomain are compatible, callers may supply the input itself as the output buffer. An implementation must therefore also be safe when out is x or out is y. If the underlying algorithm cannot naturally handle this aliasing, preserve the input through its vector space before overwriting out:

def _ieval(self, x, out):
    if out is x:
        source = self.domain.empty()
        self.domain.copyto(source, x)
        x = source
    evaluate_into(x, out)
    return out

This example uses only representation-independent storage operations. A specialized NumPy operator may instead use array assignment, while an NGSolve operator may use the corresponding NGSolve storage interface.

Cached linearization data must remain valid if the caller subsequently reuses or modifies x or out.

Defining the class

A custom linear operator typically has the following structure:

from regpy.operators import Operator

class MyOperator(Operator):
    def __init__(self, par_1, par_2):
        # Compute the domain and codomain here if they are not parameters.
        super().__init__(
            domain=my_domain,
            codomain=my_codomain,
            linear=True,
        )

    def _eval(self, x):
        # Compute y = T(x).
        return y

    def _adjoint(self, y):
        # Compute x = T.adjoint(y) for the standard pairings.
        return x

See the Volterra example for a compact application and the traction-force microscopy example for an NGSolve-based operator.

Example

Combining the preceding snippets yields the following Fourier-transform operator:

import numpy as np
from regpy.operators import Operator
from regpy.vecsps import UniformGridFcts

class SimpleFFTOnSquare(Operator):
    def __init__(self, d):
        domain = UniformGridFcts(d, d, dtype=complex)
        # Dual grid for the FFT approximation of the continuous transform.
        cd = (-1 / (2 * domain.spacing[0]),
              1 / (2 * domain.spacing[0]),
              domain.shape[0])
        codomain = UniformGridFcts(cd, cd, dtype=complex)
        super().__init__(
            domain=domain,
            codomain=codomain,
            linear=True,
        )

    def _eval(self, x):
        return np.fft.fftshift(np.fft.fftn(np.fft.ifftshift(x), norm='ortho'))

    def _adjoint(self, y):
        return np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(y), norm='ortho'))

This is a simplified version of the Fourier-transform implementation in regpy.operators.

Nonlinear operators

Nonlinear operators also define a domain, a codomain, and an _eval method. Regularization methods for inverse problems usually also need the Fréchet derivative \(F'[x]\) and its adjoint.

RegPy deliberately couples a derivative to the evaluation at its linearization point. Calling op.linearize(x) evaluates \(F(x)\), permits _eval to cache data required at \(x\), and returns both the value and a linear operator representing \(F'[x]\). Evaluating the original operator again revokes the previous derivative, because its cached data may no longer be valid. If derivatives at several points must coexist, make independent copies of the nonlinear operator before linearizing them.

The methods for evaluation, derivative and adjoint

Implement the following methods:

• _eval(x, differentiate=False) computes \(F(x)\). If differentiate is true, it must also prepare any state needed by the next two methods.

• _derivative(h) computes \(F'[x]h\) at the most recent linearization point.

• _adjoint(y) computes \(F'[x]^*y\) at that same point.

These private methods are called through op(x), op.linearize(x), and the linear operator returned by linearize; users should not invoke them directly. The corresponding in-place methods are described in +:ref:in_place_variants

The _eval method

The differentiate flag avoids derivative-specific precomputation during an ordinary evaluation. A typical allocating implementation has the form:

def _eval(self, x, differentiate=False):
    y = compute_value(x)
    if differentiate:
        self._linearization_data = prepare_derivative_data(x, y)
    return y

Store only the data needed by _derivative and _adjoint. If a stored vector must remain independent of its input, copy it through the corresponding vector space rather than assuming that the vector provides copy().

The _derivative and _adjoint method

After linearization, _derivative acts as the evaluation method of the linear operator \(F'[x]\), while _adjoint acts as its adjoint method. The returned derivative retains a revocable reference to the nonlinear operator and becomes invalid after that operator is evaluated again.

What happens when you linearize

When RegPy calls op.linearize(x), the following steps occur:

1. The operator is evaluated with differentiate=True.

2. The operator stores any intermediate quantities needed to apply the derivative and its adjoint.

3. linearize returns the value and an regpy.operators.Operator representing the derivative.

y, derivative = my_op.linearize(x)

For operators whose codomain vectors are too expensive to construct, linearize(x, return_adjoint_eval=True) returns \(F'[x]^*F(x)\) in place of \(F(x)\). The second return value is still the derivative.

Thus a typical implementation would look like this:

from regpy.operators import Operator

class MyOperator(Operator):
    def __init__(self, par_1, par_2):
        # Compute the domain and codomain here if they are not parameters.
        super().__init__(
            domain=my_domain,
            codomain=my_codomain,
            linear=False,
        )

    def _eval(self, x, differentiate=False):
        y = compute_value(x)
        if differentiate:
            self._linearization_data = prepare_derivative_data(x, y)
        return y

    def _derivative(self, h):
        return apply_derivative(self._linearization_data, h)

    def _adjoint(self, y):
        return apply_adjoint(self._linearization_data, y)

Example

As an example, consider the phase-retrieval observation operator obtained by composing a Fourier transform with the pointwise squared modulus, \(x\mapsto |\mathcal{F}(x)|^2\). Its domain is a centered complex uniform grid, and its codomain is the corresponding real Fourier grid.

def __init__(self, domain):
    # Dual grid for the FFT approximation of the continuous transform.
    cd = (-1 / (2 * domain.spacing[0]),
          1 / (2 * domain.spacing[0]),
          domain.shape[0])
    codomain = UniformGridFcts(cd, cd)
    super().__init__(
        domain=domain,
        codomain=codomain,
        linear=False,
    )

The derivative of \(Sq\colon x\mapsto |x|^2\) at \(f\) is \(h\mapsto 2\operatorname{Re}(\overline{f}h)\). Since the Fourier transform is linear, the chain rule gives

\[(Sq\circ\mathcal{F})'[f]h=Sq'[\mathcal{F}(f)]\mathcal{F}'[f]h = Sq'[\mathcal{F}(f)]\mathcal{F}(h).\]

Moreover, for the adjoint we obtain

\[(Sq\circ\mathcal{F})'[f]^\ast y= \mathcal{F}^\ast Sq'[\mathcal{F}(f)]^\ast y.\]

Both the derivative and its adjoint require \(\mathcal{F}(f)\). The evaluation method therefore stores this value only when called with differentiate=True:

def _eval(self, x, differentiate=False):
    y = np.fft.fftshift(np.fft.fftn(np.fft.ifftshift(x), norm='ortho'))
    if differentiate:
        self._factor = y
    return y.real**2 + y.imag**2

The _derivative method applies this derivative to a direction \(h\). Because \(\mathcal{F}(f)\) is already cached, the implementation is:

def _derivative(self, h):
    return 2*(self._factor.conj() * np.fft.fftshift(np.fft.fftn(np.fft.ifftshift(h), norm='ortho'))).real

The adjoint maps \(y\) to \(\mathcal{F}^\ast(2\mathcal{F}(f)y)\) and reuses the same cached factor:

def _adjoint(self, y):
    return np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(2*self._factor * y), norm='ortho'))

Together, these methods define the phase-retrieval observation operator:

import numpy as np
from regpy.vecsps import UniformGridFcts
from regpy.operators import Operator

class Observation(Operator):

    def __init__(self, domain):
        # Dual grid for the FFT approximation of the continuous transform.
        cd = (-1 / (2 * domain.spacing[0]),
              1 / (2 * domain.spacing[0]),
              domain.shape[0])
        codomain = UniformGridFcts(cd, cd)
        super().__init__(
            domain=domain,
            codomain=codomain,
            linear=False,
        )

    def _eval(self, x, differentiate=False):
        y = np.fft.fftshift(np.fft.fftn(np.fft.ifftshift(x), norm='ortho'))
        if differentiate:
            self._factor = y
        return y.real**2 + y.imag**2

    def _derivative(self, h):
        return 2*(self._factor.conj() * np.fft.fftshift(np.fft.fftn(np.fft.ifftshift(h), norm='ortho'))).real

    def _adjoint(self, y):
        return np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(2*self._factor * y), norm='ortho'))

In practice, this operator need not be implemented from scratch. RegPy can compose the existing Fourier-transform and squared-modulus operators:

from regpy.operators import SquaredModulus, FourierTransform
from regpy.vecsps import UniformGridFcts

domain = UniformGridFcts((-1, 1, 100), (-1, 1, 100), dtype=complex)

ft = FourierTransform(domain, centered=True)
sqm = SquaredModulus(ft.codomain)

observe = sqm * ft

Further examples

See the Volterra example for another nonlinear operator.

NGSolve operators

RegPy integrates with NGSolve through regpy.operators.ngsolve, regpy.hilbert.ngsolve, and regpy.functionals.ngsolve.

As for every RegPy operator, adjoints are implemented with respect to the standard vector-space pairings, not a separately chosen Hilbert-space inner product.

Use regpy.operators.ngsolve.NgsOperator as the base class for custom NGSolve operators.

Caution

The NGSolve interface is still evolving and may change.

For coefficient identification in second-order elliptic PDEs, RegPy provides regpy.operators.ngsolve.SecondOrderEllipticCoefficientPDE. Subclasses specify the required bilinear and linear forms. See the diffusion-coefficient example.

Combined adjoint and derivative

Some operators admit a more efficient implementation of \(F'[x]^*F'[x]\) than separate derivative and adjoint applications. This is particularly useful when codomain vectors would be too large to store. Such an operator may implement _adjoint_derivative directly:

def _adjoint_derivative(self, h):
    return apply_normal_operator(self._linearization_data, h)

Warning

Support for combined adjoint-derivative evaluations is specialized. Check that the intended solver uses this interface before relying on it.