1 Introduction↩︎

In this paper we are concerned with inverse problems for which the forward operator \(K\) has a significant null space. More precisely, we assume that \(K: L^2(\Omega) \rightarrow L^2(\partial\Omega)\) is a linear operator and consider the weighted problem \[\label{def:variational95form} \min_{f \in BV(\Omega)} \left\{ \frac{1}{2} \| Kf - d \|_{L^2(\partial\Omega)}^2 + \alpha \underbrace{\int_{\Omega} w(y) \, |Df| (y)}_{= \operatorname{TV}_w(f)} + \beta \int_{\partial\Omega} w_\partial (y) |f(y)| \, ds(y) \right\}.\tag{1}\] Here, \(d \in L^2(\partial\Omega)\) is the measured data, \(\Omega \subset \mathbb{R}^2\) is a bounded Lipschitz domain, \(\alpha > 0\) is a regularization parameter, \(Df\) is the vector-valued Radon measure representation of the "gradient" of \(f\), and \(w\) is the weight function \[\label{def:w} w(y)=\sup_{\boldsymbol{d},\; |\boldsymbol{d}|_2 =1} \left\| K\!\left( \nabla G(\cdot;y)\cdot \boldsymbol{d} \right) \right\|_{L^p(\partial\Omega)},\tag{2}\] where \(G\) is a suitable Green’s function which we define below. The presence of the boundary term, i.e., the term involving the parameter \(\beta\), is motivated by the following observation: Without this term, it becomes "cheap" to obtain a good data fit with a source which is (approximately) constant close to \(\partial \Omega\) since the total variation of a constant function is zero. Hence, a boundary-bias occurs. We will return to this issue below and explain how it is linked to the choice of boundary conditions to use in the definition of the Green’s function \(G\).

Operators that map interior distributions to boundary observations, i.e., \(K \colon L^2(\Omega) \to L^2(\partial \Omega)\), arise naturally in inverse source problems for elliptic partial differential equations, see [1], [2] and the references therein. Concrete examples include source identification in groundwater modelling, the inverse problem of electrocardiography [3], and certain regimes of electrical impedance tomography [4]. A defining feature of such operators is that the dimension of the null space is typically infinite: many distinct interior sources produce identical boundary data. Stable reconstruction therefore requires regularization, and the regularization functional plays a particularly active role since it must determine the component of the solution that lies in \(\mathcal{N}(K)\).

When \(K\) has a substantial null space, the conventional choice \(w \equiv 1\) in 1 is known to introduce systematic artifacts in the reconstructed solution. Since each unit of total variation produces a far larger reduction of the data-fidelity term in regions where \(K\) is highly sensitive than in regions where \(K\) is nearly blind, an unweighted TV penalty effectively becomes cheaper, relative to the fidelity term, in the sensitive regions. As a consequence, recovered features tend to be displaced toward those regions, and sources of equal magnitude may be reconstructed with different intensities - if at all - depending on their location relative to \(\partial \Omega\). These observations motivate spatially varying weights \(w\) that rebalance the TV penalty according to the non-uniform sensitivity of \(K\).

Spatially weighted TV functionals have been considered by several authors and in several contexts. Anisotropic and direction-aware TV functionals are studied in [5], [6], and structure-guided variants tailored to multi-channel and multi-modal data appear in [7], [8]. Specifically for inverse problems with significant null spaces, weighted Tikhonov and sparsity formulations have been investigated previously, and the anisotropic (i.e., \(|Df|_1\)) weighted TV counterpart to problem 1 was introduced and analyzed in [9]. The present paper builds on that work in two main respects: (i) it replaces the axis-aligned weights in the anisotropic setting with taking the supremum over all possible directions, and (ii) it introduces and analyzes the directionally weighted refinement \(\operatorname{TV}_{\widetilde{w}}\), in which the worst-case supremum over directions is replaced by the actual jump direction of the candidate solution, see 3 and 4 below.

The form of the weight 2 can be motivated as follows. Using the representation of \(K\) in terms of \(G\), the action of \(K\) on a localized perturbation of \(f\) near a point \(y\), in direction \(\boldsymbol{d}\), is to leading order given by \(K(\nabla G(\cdot;y)\cdot \boldsymbol{d})\). The weight \(w(y)\) thus measures the maximal sensitivity of the data to a unit perturbation of the gradient of \(f\) at \(y\), with the supremum taken over all admissible directions. The role of \(w\) in 1 is thus to remove spatial bias: where \(K\) is weakly sensitive, \(w(y)\) is small and the TV penalty becomes inexpensive, so the minimizer is allowed to develop jumps and other structure there. Conversely, in regions of high sensitivity each unit of total variation already produces a large change in \(Kf\), so a larger value of \(w(y)\) is needed to prevent the data-fidelity term from biasing reconstructions toward those regions. In this way, the weight rebalances the cost of total variation against the heterogeneous sensitivity of \(K\) across \(\Omega\), allowing the reconstruction to follow the data wherever the data are informative and the prior wherever they are not.

Our investigation will also suggest the following alternative to 1 , involving the polar decomposition \(Df = {\boldsymbol{\nu}}_f |Df|\), \(|{\boldsymbol{\nu}}_f|_2 = 1\), \[\label{def:variational95form95pre} \min_{f \in BV(\Omega)} \left\{ \frac{1}{2} \| Kf - d \|_{L^2(\partial)}^2 + \alpha \underbrace{\int_{\Omega} \widetilde{w}(y,f) \, |Df|(y)}_{= \operatorname{TV}_{\widetilde{w}}(f)} + \beta \int_{\partial\Omega} w_\partial (y) |f(y)| \, ds(y)\right\},\tag{3}\] where the weight function in this case reads \[\label{def:tw} \widetilde{w}(y,f) = \left\| K\!\left( \nabla G(\cdot;y)\cdot {\boldsymbol{\nu}}_f(y) \right) \right\|_{L^p(\partial\Omega)}.\tag{4}\] Below we will simply write \({\boldsymbol{\nu}}\), instead of \({\boldsymbol{\nu}}_f\), whenever it is clear from the context with which function \(f\), \({\boldsymbol{\nu}}\) is associated. The crucial difference between 2 and 4 is that the supremum over \(\boldsymbol{d}\) in 2 is replaced by the actual jump direction \({\boldsymbol{\nu}}(y)\) of the candidate \(f\). While 2 treats every direction equally and thereby penalizes the worst-case sensitivity at \(y\), the formulation in 4 is sensitive to the direction of the level curves: variations of \(f\) in different directions are penalized differently, in accordance with how strongly \(K\) responds to perturbations in those directions. Setting \(\phi(y, {\boldsymbol{\nu}}) := \|K(\nabla G(\cdot;y) \cdot {\boldsymbol{\nu}})\|_{L^p(\partial\Omega)}\), the map \({\boldsymbol{\nu}}\mapsto \phi(y,{\boldsymbol{\nu}})\) is the composition of a linear map with the \(L^p\) norm and is therefore positively \(1\)-homogeneous and sublinear on \(\mathbb{R}^2\) for every \(y\). Consequently, \[\operatorname{TV}_{\widetilde{w}}(f) \;=\; \int_\Omega \phi(y, {\boldsymbol{\nu}}(y)) \, d|Df|(y)\] fits into the framework of anisotropic, position-dependent total variation functionals (i.e., Finsler TV [10]) \[\label{def:phi95TV} f \mapsto \int_\Omega \phi(y, dDf(y)),\tag{5}\] in which \(\phi: \overline{\Omega} \times \mathbb{R}^N \to [0,\infty)\) is positively \(1\)-homogeneous in its second argument; despite the appearance of \({\boldsymbol{\nu}}\) in the integrand, \(\operatorname{TV}_{\widetilde{w}}(f)\) is just the polar-decomposition representation of 5 for our specific choice of \(\phi\). Functionals of this type are convex and lower semicontinuous on \(BV(\Omega)\) under mild assumptions on \(\phi\), and have been studied since the 1990s in the context of relaxation and integral representation [11]–[13], as well as in connection with weighted-anisotropic TV regularization for image reconstruction [14]. The class 5 continues to attract attention: the gradient flow associated with general convex linear-growth integrands \(f(x,\mathbb{A}u)\) and the corresponding BV-relaxation are studied in [15], and a general lower semicontinuity and existence theory for anisotropic TV functionals with position-dependent integrand is developed in [16]. The closely related directional and structure-guided TV functionals appear in [6], [7].

Finally, by construction one has the pointwise inequality \(\operatorname{TV}_{\widetilde{w}}(f) \le \operatorname{TV}_w(f)\), with equality at points where the jump direction \({\boldsymbol{\nu}}(y)\) realizes the supremum in 2 .

In addition to 1 and 3 , we also consider their basis-pursuit counterparts, in which the data fidelity term is replaced by a constraint: \[\min_{f \in BV(\Omega)} \operatorname{TV}_w(f) \quad \text{subject to} \quad K f = d,\] and analogously for \(\operatorname{TV}_{\widetilde{w}}\).

1.0.0.1 Contributions and outline.

The purpose of this paper is to (i) derive the weights \(w\) and \(\widetilde{w}\) from a sensitivity argument involving the Green’s function \(G\) associated with \(K\), (ii) establish some recoverability properties of the basis pursuit counterparts to 1 and 3 , and (iii) illustrate, by means of numerical experiments, that the weighted formulations correct several of the artifacts produced by standard isotropic TV when \(K\) has a significant null space.

The remainder of the paper is organized as follows. Section 2 provides a motivating example where unweighted total variation fails to recover an interior source. Section 3 collects the necessary background on the Green’s function setting, and the detailed motivation and analysis of \(w\) and \(\widetilde{w}\). In Section 4 we present the optimality conditions and establish some recoverability criteria for the basis-pursuit counterparts of 1 and 3 . Section 5 reports numerical experiments illustrating the theoretical results, and Section 6 contains concluding remarks.

2 Motivating example↩︎

Here we present a numerical example showing that standard unweighted isotropic TV fails to yield adequate results for identifying an internal source from boundary measurements. Throughout this section we let \(\Omega = (0,1)^2\) and take \(K\) to be the operator that maps a source \(f \in L^2(\Omega)\) to the boundary trace \(u|_{\partial\Omega}\) of the solution \(u \in H^1(\Omega)\) of the screened-Poisson problem \[\label{eq:screenedPoisson} -\Delta u + u = f \quad \text{in } \Omega, \qquad u = 0 \quad \text{on } \partial \Omega,\tag{6}\] for which the analysis below applies. Boundary data \(d = K f^\dagger\) are generated synthetically from a known/true piecewise constant source \(f^\dagger\), and we then solve 3 with \(\widetilde{w}\equiv 1\). The true source is the indicator function of a disk centered at \((0.5, 0.6)\) with radius \(0.3\), see Figure 1 (a).

The reconstructions obtained with unweighted TV are shown for \(\beta = 0\) (i.e., without boundary penalty) in Figure 1 (b) and for \(\beta = \alpha\) (i.e., with boundary penalty) in Figure 1 (c). Two effects are immediately visible. First, the recovered support is pulled toward the boundary \(\partial \Omega\): the reconstructed mass is concentrated where it produces the largest boundary signal per unit total variation, rather than at the true location of the source. This is particularly outspoken when \(\beta = 0\). Second, the magnitude of the reconstruction is markedly smaller than that of \(f^\dagger\). Both effects originate from the spatial inhomogeneity of \(K\): in the unweighted formulation, every unit of \(|Df|\) is charged the same price regardless of where it is placed, even though a unit of \(|Df|\) near \(\partial\Omega\) explains far more of \(d\) than a unit of \(|Df|\) deep in the interior. The minimizer therefore prefers configurations near the boundary, and accepts a global reduction in amplitude in order to lower the regularization cost.

This example sets the stage for the rest of the paper. It shows that, when the forward operator has a significant null space, the choice of regularization functional must account for the spatial sensitivity of \(|Df|\), which serves as the main motivation for introducing the weight functions \(w\) and \(\widetilde{w}\).

3 Preliminaries↩︎

We will now derive how the source \(f\) can be represented in terms of the Green’s function of the Laplace operator and the gradient \(\nabla f\). These standard results then motivate the definition of our weights and can be used to establish some basic inequalities.

One interesting point is that the choice of boundary conditions for the Green’s function influences the weights and the choice could - and perhaps should - be based on whether one wants to penalize jumps at the boundary or not, i.e., whether \(\beta = 0\) or \(\beta > 0\) in 1 and 3 . We will return to this after introducing the framework.

Neumann representation↩︎

Let \(G\) denote the Green’s function associated with a standard self-adjoint second order elliptic operator subject to homogeneous Neumann boundary conditions: \[\label{def:Greens95function} \begin{align} &-\Delta_y G(y;x)=\delta_x(y)-\frac{1}{|\Omega|} \quad y \in \Omega, \\ &\nabla_y G(y;x) \cdot \boldsymbol{n}(y) = 0 \quad y \in \partial \Omega, \end{align}\tag{7}\] and note that we may express the source \(f\) in the form, using integration by parts, \[\begin{align} \nonumber f(x) &= \int_\Omega \delta_x(y) f(y) \, dy \\ \nonumber &=\int_{\Omega} \nabla G(y;x) \cdot Df(y) + \frac{1}{|\Omega|} \int_\Omega f(y) \, dy \\ \nonumber &=\int_{\Omega} \nabla G(y;x) \cdot {\boldsymbol{\nu}}(y) \, |Df|(y) + \frac{1}{|\Omega|} \int_\Omega f(y) \, dy\\ \label{eq:source95representation} &=\int_{\Omega} \nabla G(x;y) \cdot {\boldsymbol{\nu}}(y) \, |Df|(y) + \frac{1}{|\Omega|} \int_\Omega f(y) \, dy, \end{align}\tag{8}\] where we have employed the polar decomposition \(Df = {\boldsymbol{\nu}}|Df|\), \(| {\boldsymbol{\nu}}|_2 = 1\), and the symmetry of the Green’s function \(G\). We also observe that this representation is independent of the free constant of the solution of 7 .

Let us, for the sake of simplicity, assume that \(K1=0\), i.e., that the null space of \(K\) contains the constant functions. If not, we could consider the composite function \(Q = P \circ K\), where \(P\) projects any \(y \in R(K)\) onto the orthogonal complement of the image under \(K\) of the constants and then define the weights in 2 and 4 using \(Q\) rather than \(K\). See [9] for details.

From 8 we find that, keeping in mind that \(K\) is a linear operator, \[\begin{align} Kf &= K \int_{\Omega} \nabla G(\cdot;y) \cdot {\boldsymbol{\nu}}(y) \, |Df|(y) \\[0.5em] &= \int_{\Omega} K\!\left( \nabla G(\cdot;y)\cdot {\boldsymbol{\nu}}(y) \right) \, |Df|(y) , \end{align}\] and we get the following bound for the \(L^p\)-norm of the image of \(f\) under \(K\): \[\begin{align} \nonumber \|Kf\|_{L^p(\partial\Omega)} &= \left(\int_{\partial\Omega} \left| \int_\Omega K\!\left( \nabla G(\cdot;y)\cdot {\boldsymbol{\nu}}(y) \right) \, |Df|(y) \right|^p ds(z) \right)^{1/p} \\[0.5em] \nonumber &\le \int_{\Omega} \left( \int_{\partial\Omega} \left| K\!\left( \nabla G(\cdot;y)\cdot {\boldsymbol{\nu}}(y) \right) \right|^p ds(z) \right)^{1/p} \, |Df|(y) \\[0.5em] \nonumber &= \int_{\Omega} \left\| K\!\left( \nabla G(\cdot;y)\cdot {\boldsymbol{\nu}}(y) \right) \right\|_{L^p(\partial\Omega)} \, |Df|(y) \\[0.5em] \label{eq:upper95bound95pre} &= \operatorname{TV}_{\widetilde{w}}(f), \end{align}\tag{9}\] provided that \(\widetilde{w}\) is defined according to 4 .

It follows from the definitions 2 and 4 that \[\widetilde{w}(y,f) \leq w(y) \quad \forall y \in \Omega, \, \forall f \in BV(\Omega),\] and hence that \[\operatorname{TV}_{\widetilde{w}}(f) \leq \operatorname{TV}_w(f) \quad \forall f \in BV(\Omega).\] We can therefore in view of 9 conclude that also \[\begin{align} \label{eq:upper95bound} \|Kf\|_{L^p(\partial\Omega)} &\leq \operatorname{TV}_w(f) \quad \forall f \in BV(\Omega). \end{align}\tag{10}\]

Dirichlet representation↩︎

If we instead consider the Green’s function for the Laplace operator with homogeneous Dirichlet boundary conditions, we get, in contrast to 8 , an alternative representation of the source \(f\), namely \[\label{label:source95rep95dirichlet} f(x) = \int_\Omega \nabla G(x;y) \cdot {\boldsymbol{\nu}}(y)|Df|(y) - \int_{\partial\Omega} \partial_{\boldsymbol{n}}G(x;y) \;f(y) \, ds(y).\tag{11}\]

Assuming \(f\) to have zero trace, we get exactly the same upper bounds as derived in 9 and 10 . If not, we get, from an argument analogous to 9 , that \[\|Kf\|_{L^p(\partial\Omega)} \leq \operatorname{TV}_{\widetilde{w}}(f) + \int_{\partial\Omega} w_\partial |f| \, ds \leq \operatorname{TV}_w(f) + \int_{\partial\Omega} w_\partial |f| \, ds,\] where \[w_\partial(y) = \left(\int_{\partial\Omega} |K (\partial_{\boldsymbol{n}} G(\cdot;y))|^p ds(z)\right)^{1/p}.\]

4 Analysis↩︎

In this section we will first consider some instructive cases where we can guarantee exact recovery of a given source, before we view the problems through the lens of more classical optimality conditions.

4.1 Exact recovery↩︎

We will now present an exact recovery result for the zero regularization counterpart to 3 :

The following quantity plays an important role in our analysis of 1 , \[\label{def:optimal95direction} \boldsymbol{q}(y) \in \mathop{\mathrm{arg\,max}}_{\boldsymbol{d},\;\|\boldsymbol{d}\|=1} \left\| K\!\left(\nabla G(\cdot;y)\cdot \boldsymbol{d} \right) \right\|_{L^1(\partial \Omega)} \quad y \in \Omega .\tag{12}\] We remark that \(\boldsymbol{q}(y)\) is not unique and that Theorem 3 below holds for any \(\boldsymbol{q}(y)\) satisfying 12 , provided that also ?? ?? hold. Roughly speaking, \(\boldsymbol{q}(y)\) is a dominant direction of \(\nabla G(\cdot;y)\) under the image of \(K\). Using this concept, our exact recovery result for the basis pursuit problem associated with 1 reads:

Example 1↩︎

Let us consider a 1D problem with \(\Omega=(0,1), \;p = 1, \;\beta = 0\) and assume that \(f^*\) is strictly increasing on \(\Omega\). Then \(|Df^*|(y) = (Df^*)(y)\) for all \(y \in \Omega\). Consequently, \[{\boldsymbol{\nu}}^*(y)= 1 \in \mathop{\mathrm{arg\,max}}_{d,\;\|d\|=1} \left\| K\!\left(G'(\cdot;y)\cdot d \right) \right\|_{L^1(\partial \Omega)} \quad y \in \Omega ,\] and it follows that ?? holds with \(q(y)=1\), \(y\in \Omega\). We conclude that ?? holds for strictly increasing functions provided \(K\) and \(G\) are such that ?? or ?? is satisfied with \({\boldsymbol{\nu}}^*(y)=1\). The same argument can, of course, be made for strictly decreasing functions.

Example 2↩︎

We consider a true source in 1D with a single jump: \[f^*(x) = \left\{ \begin{array}{cc} r-1, & x<r, \\ r, & x>r, \end{array}\right. . \label{eq:heaviside}\tag{13}\] where \(r \in (0,1)\) is fixed. Also, \(\Omega=(0,1), \;p = 1\) and \(\beta = 0\). Then, cf. 10 , \[\begin{align} \nonumber \|Kf^*\|_{L^1(\partial\Omega)} &= \int_{\partial\Omega} \left| \int_\Omega K\!\left( \nabla G(\cdot;y)\cdot {\boldsymbol{\nu}}^*(y) \right) \, |Df^*|(y) \right| dz \\[0.5em] &= \int_{\partial\Omega} \left| \int_\Omega K\!\left(G'(\cdot;y)\cdot \nu^*(y) \right) \, d \delta_r (y) \right| dz \\ &= \int_{\partial\Omega} \left| K\!\left(G'(\cdot;r)\cdot 1 \right) \right| dz \\ &=w(r). \end{align}\] Moreover, \[\begin{align} TV_w(f^*) &= \int_{\Omega} w(y) \, |Df^*| (y) \\ &= \int_{\Omega} w(y) \, d \delta_r (y) \\ &= w(r). \end{align}\] Consequently, if \(Kg=Kf^*\), then \[\begin{align} TV_w(f^*)=w(r)=\|Kf^*\|_{L^1(\partial \Omega)}=\|Kg\|_{L^1(\partial \Omega)} \leq TV_w(g), \end{align}\] where we have used 10 , and it follows that \(f^*\) obeys ?? . In this example \({\boldsymbol{\nu}}^*(y)=0\) for \(y \neq r\) and \({\boldsymbol{\nu}}^*(r)=1\) and therefore ?? , ?? and ?? all hold.

4.2 Optimality conditions↩︎

The optimality conditions associated with the noise-free limit version of 3 are presented in the proposition below. To this end, we define the polar \(F^\circ\) of a convex, lower semi-continuous and 1-homogeneous function \(F\) by \[F^\circ(y, \Psi(y)) = \sup_{\theta \neq 0} \left\{\frac{\theta \cdot \Psi(y)}{F(y,\theta)}: \theta \in \mathbb{R}^2 \right\}\]

 1.  ***Interior:**
     $\;\;\;\widetilde{w}^\circ(y,\boldsymbol{z}(y)) \leq 1, \quad y \in \Omega$,*

 2.  ***Boundary:**
     $|\boldsymbol{z}(y) \cdot \boldsymbol{n}| \leq \beta w_\partial(y), \quad y \in \partial\Omega .$*

  1.  ***Interior alignment:**
      $$\boldsymbol{z} \cdot {\boldsymbol{\nu}}_{f^*} = \widetilde{w}(y, {\boldsymbol{\nu}}_{f^*}), \quad |D f^*|-a.e.$$*

  2.  ***Boundary alignment:**
      $(\boldsymbol{z} \cdot \boldsymbol{n})f^* = - \beta w_\partial|f^*| \quad \textrm{on} \;\partial\Omega.$*

Note that if \(\beta = 0\) the dual constraint for the boundary implies that \(\boldsymbol{z} \cdot \boldsymbol{n} = 0\) along \(\partial \Omega\).

From these considerations, we get the simplified optimality conditions

This can be interpreted as follows: The TV atom \(\chi_R\) associated with a simple set \(R\) (cf. [19], [20]) is generated by a continuum of infinitesimal oriented boundary elements \((y_0,\theta_0)\) of its perimeter. Under the weight \(\widetilde{w}\) these elements are treated without spatial bias in the sense that every \((y_0,\theta_0)\) is equally dual-feasible through its own unit-length canonical certificate/Lagrange multiplier \(\lambda_{(y_0,\theta_0)}\).

Admittedly, to certify a full region, i.e., to recover \(f^* = \chi_{R}\), we would need a saturating dual certificate in a \(|Df^*|\)-a.e. sense. The existence of such a certificate is far from trivial and relies on the domain \(\Omega\), the shape \(R\) and the forward operator \(K\). In the unweighted setting, exact recovery of this type has been established for certain shapes under specific forward operators, both for isotropic TV [21] and in a anisotropic TV setting [22].

Nevertheless, we will attempt to further highlight why the weights are necessary to certify extended shapes at different depths throughout the domain \(\Omega\). To this end, assume the true source to be the characteristic function \(\chi_E\), where \(E \subset \subset \Omega\) has a smooth boundary \(\partial E\).

Note that on characteristic functions the weighted total variation reduces to the weighted perimeter, \[\operatorname{TV}_{\widetilde{w}}(\chi_E)=\int_{\partial E} \widetilde{w}(y,\boldsymbol{n}_E)ds.\] Consequently, if \(\chi_E\) solves ?? , it is a stationary point of the Lagrangian \[\mathcal{L}(E,\lambda)=\int_{\partial E}\widetilde{w}(y,\boldsymbol{n}_E)\,ds+\langle\lambda,K\chi_E-d\rangle,\] where \(\lambda\) is an associated Lagrange multiplier.

To derive the associated Euler–Lagrange equation, consider a smooth normal deformation \[\partial E_t = \{\, y + t\,\varphi(y)\,\boldsymbol{n}_E(y) : y \in \partial E \,\},\] with \(\varphi\) an arbitrary smooth scalar velocity. The stationarity condition of the Lagrangian reads \[\frac{d}{dt}\, \mathcal{L}(E_t,\lambda)\bigg|_{t=0} = 0 \quad \forall\, \varphi .\] By the theory of shape derivatives [23], the first variation of the constraint term is \[\frac{d}{dt}\, \langle \lambda,\, K\chi_{E_t}\rangle \bigg|_{t=0} = \int_{\partial E} K^*\lambda \;\varphi \; ds,\] while the first variation of the weighted perimeter (in which the normal \(\boldsymbol{n}_{E_t}\) also varies with \(t\)) is \[\frac{d}{dt} \int_{\partial E_t} \widetilde{w}(y, \boldsymbol{n}_{E_t})\, ds \bigg|_{t=0} = \int_{\partial E} \kappa_{\widetilde{w}}\; \varphi \; ds,\] where \(\kappa_{\widetilde{w}}\) is the weighted anisotropic curvature [24]. Stationarity for all \(\varphi\) therefore gives \[\int_{\partial E} \big(\kappa_{\widetilde{w}} + K^*\lambda\big)\,\varphi \; ds = 0 \quad \forall\, \varphi,\] and hence \[\label{eq:euler95lagrange95shape} \kappa_{\widetilde{w}} + K^*\lambda = 0 \quad \text{on } \partial E .\tag{14}\] We emphasize that 14 is a necessary (stationarity) condition that must hold simultaneously at every point of \(\partial E\); it does not by itself guarantee that \(\chi_E\) is a global minimizer.

For simplicity, assume now that \(\widetilde{w}(y,{\boldsymbol{\nu}}) = w(y)\), i.e.that the weight is isotropic and that \(E = B_r(x_0)\), i.e., a ball of radius \(r\) centered at \(x_0 \in \Omega\). The weighted curvature then simplifies to \[\kappa_{\widetilde{w}} = \frac{w}{r} + \nabla w \cdot \boldsymbol{n},\] so that 14 becomes \[\label{eq:curvature95cond} K^*\lambda = -\frac{w}{r} - \nabla w \cdot \boldsymbol{n} \quad \text{on } \partial B_r(x_0) .\tag{15}\]

Let us now discuss the realism of satisfying 15 for different choices of weight functions \(w\). To this end, let us define \(\rho(y) := \operatorname{dist}(y,\partial\Omega)\) as the distance from \(y\) to the boundary, which we refer to as the depth of \(y\). We also assume, which is the case for many inverse problems, that \(K^*\lambda\) decays exponentially into the interior: \[|K^*\lambda(y)| \;\lesssim\; C(\rho(y))\, e^{-\tau \rho(y)},\] where \(C(\cdot)\) has at most polynomial dependence on \(\rho\).

Note that 15 must hold at every point of \(\partial B_r\) for the same multiplier \(\lambda\). At a single point, the equation can always be satisfied; the difficulty is simultaneity along the curve \(\partial B_r\). With unweighted TV, i.e., \(w \equiv 1\), the right hand side becomes the constant \(-\frac{1}{r}\), which cannot match the (assumed) exponentially decay of \(K^*\lambda\) for curves spanning over multiple depths. A weight varying only polynomially in \(\rho\) also fails for the same reason. Hence, to allow an approximate recovery of interfaces that span a range of depths, it appears that the weight function must have the same exponential dependence on \(\rho\) as \(K^*\lambda\). This is achieved in the present paper by defining \(w\) in terms of images under \(K\) of the gradient of appropriate Green’s functions.

5 Numerical results↩︎

Discretization↩︎

All experiments are carried out on the unit square \(\Omega = (0,1) \times (0,1)\), discretized by a uniform grid with spacing \(h = 1/128\). The forward operator \(K: L^2(\Omega) \to L^2(\partial\Omega)\) is the map \(f \mapsto u|_{\partial\Omega}\), where \(u\) solves the screened Poisson equation 6 . The PDE is discretized with standard piecewise-linear finite elements on the triangulation associated with the grid, and the resulting discretized forward operator is denoted by \(\mathsf{K} = \mathsf{T}\mathsf{A}^{-1} \in \mathbb{R}^{M \times N}\), where \(\mathsf{A}\) and \(\mathsf{T}\) are the matrices associated with the screened Poisson equation and the trace operator, respectively.

To describe how the weights are generated, let \(\mathsf K = \mathsf U\mathsf S\mathsf V^\top\) be the (truncated) singular value decomposition of the forward matrix, with singular values \(s_1,\dots,s_r\) and right singular vectors \(\boldsymbol{v}_1,\dots,\boldsymbol{v}_r\in\mathbb{R}^N\). Furthermore, we denote by \(\mathsf{L}\) the discretization of the Laplace operator \(-\Delta\), used for the representation of the Green’s functions, cf. 11 . The directional weight at a point \(y\) is the magnitude of the dipole sensitivity, \[\label{eq:sensitivity} s_{\boldsymbol{\nu}}(y) = \bigl\| \mathsf K\,\nabla_{\boldsymbol{\nu}} G(\cdot\,;y) \bigr\|_2, \qquad G(\cdot\,;y) = \mathsf L^{-1}\delta_y ,\tag{16}\] i.e., the discrete counterpart of \(\widetilde{w}(y,\boldsymbol{\nu}) = \|K(\nabla G(\cdot\,;y)\cdot\boldsymbol{\nu})\|_{L^2(\partial \Omega)}\).

Evaluating 16 directly would require one dipole solve per cell, i.e.\(N\) solves. Index the grid cells by \(e\), with \(y_e\) the corresponding point and \(\delta_e := \delta_{y_e}\) the discrete source at that cell. The directional derivative \(\nabla_{\boldsymbol{\nu}}\) acts on the source location \(y\), while \(\mathsf{L}^{-1}\) acts on the field variable; thus linearity of \(\mathsf{L}^{-1}\) gives \(\nabla_{\boldsymbol{\nu}}\mathsf{L}^{-1}\delta_e = \mathsf{L}^{-1}\nabla_{\boldsymbol{\nu}}\delta_e\). Using in addition the symmetry \(\mathsf L = \mathsf L^\top\) to write \(\boldsymbol{v}_k^\top\mathsf{L}^{-1} = \boldsymbol{z}_k^\top\) with \(\boldsymbol{z}_k = \mathsf{L}^{-1}\boldsymbol{v}_k\) gives \[\boldsymbol{v}_k^\top \mathsf{L}^{-1}\nabla_{\boldsymbol{\nu}}\delta_e = (\mathsf{L}^{-1}\boldsymbol{v}_k)^\top\nabla_{\boldsymbol{\nu}}\delta_e = \boldsymbol{z}_k^\top \nabla_{\boldsymbol{\nu}}\delta_e = (\nabla_{\boldsymbol{\nu}}^\top\boldsymbol{z}_k)_e = -(\nabla_{\boldsymbol{\nu}}\boldsymbol{z}_k)_e,\] where the last equality follows from integration by parts. Expanding 16 in the singular basis of \(\mathsf{K}\) and using the expression above yields \[\label{eq:reciprocity} \bigl\| \mathsf K\,\nabla_{\boldsymbol{\nu}} G(\cdot\,;y_e)\bigr\|_2^2 = \sum_{k=1}^{r} s_k^2\, \bigl(\boldsymbol{v}_k^\top \mathsf L^{-1}\,\nabla_{\boldsymbol{\nu}}\delta_e\bigr)^2 = \sum_{k=1}^{r} s_k^2\,\bigl(\nabla_{\boldsymbol{\nu}}\,\boldsymbol{z}_k\bigr)_e^2.\tag{17}\] Thus pushing each detectable mode/singular vector once through \(\mathsf L^{-1}\) and reading off its discrete gradient recovers the sensitivity at every cell and in every direction simultaneously, replacing \(N\) per-cell solves by \(r\) mode solves (\(r\ll N\)).

Stacking the gradients of the smoothed modes column-wise gives \[{\mathsf W}_x = \mathsf D_x\,[\,\boldsymbol{z}_1\;\cdots\;\boldsymbol{z}_r\,]\,\mathsf S_r, \qquad {\mathsf W}_y = \mathsf D_y\,[\,\boldsymbol{z}_1\;\cdots\;\boldsymbol{z}_r\,]\,\mathsf S_r, \qquad \mathsf S_r = \operatorname{diag}(s_1,\dots,s_r),\] with \(\mathsf D_x, \mathsf D_y\) the discrete partial-derivative operators. For each evaluation cell \(e\) (paired with one \(x\)- and one \(y\)-edge) the corresponding rows form the local weight \[{\mathsf W}^{(e)} = \begin{bmatrix} ({\mathsf W}_x)_{e,:} \\[2pt] ({\mathsf W}_y)_{e,:} \end{bmatrix}^{\!\top} \in \mathbb{R}^{r\times 2},\] whose \(2\times 2\) Gram matrix is precisely the metric whose quadratic form in \(\boldsymbol{\nu}\) equals \(s_{\boldsymbol{\nu}}(y_e)^2\) in 16 .

With \((\nabla_h \boldsymbol{f})_e = \bigl((\mathsf D_x \boldsymbol{f})_e,(\mathsf D_y \boldsymbol{f})_e\bigr)^\top\), the discrete directionally weighted total variation (\(p=2\)) is the group (\(\ell_{2,1}\)) norm of the weighted gradient, \[\mathrm{TV}_{\widetilde{w}}(\boldsymbol{f}) = \sum_{e=1}^{N_{\mathrm{eval}}} \bigl\| {\mathsf W}^{(e)}(\nabla_h \boldsymbol{f})_e \bigr\|_2 .\]

From this, it is also straightforward to compute the discretized isotropic weights corresponding to 2 by setting \[\mathrm{TV}_{w}(\boldsymbol{f}) = \sum_{e=1}^{N_{\mathrm{eval}}} \|\mathsf W^{(e)}\|_2 \,\bigl\|(\nabla_h \boldsymbol{f})_e\bigr\|_2 .\]

Equation 17 shows that each mode’s contribution \((\nabla_{\boldsymbol{\nu}}\boldsymbol{z}_k)_e^2\) involves the factor \(s_k^2\). In our experiments, removing this factor - i.e. setting \(s_k = 1\) for all \(k\) in 17 - improves recovery in practise, and all numerical results below use this choice. For a more in-depth discussion of employing such a filter for weight generation, see [25]. Note that our analysis is presented in terms of an abstract forward operator \(K\). Hence, such a filter can also be incorporated in the infinite dimensional setting by considering it to be part of the action of \(K\).

Test problems↩︎

We consider two test problems, which differ in the geometry of the unknown. The first is an inverse source problem: the goal is to recover a localized source \(f^\dagger\) supported in the interior of \(\Omega\) from boundary measurements of the corresponding solution \(u\). The ground truth has compact support strictly inside \(\Omega\), so feature creation on \(\partial \Omega\) is undesirable. The second is a gravimetry-inspired problem: the unknown \(f^\dagger\) models a piecewise constant density and the discontinuity of interest is the interface between two materials of different density, which need not stay away from the boundary.

This distinction directly informs the choice of Green’s function and, with it, the form of the regularizer. As discussed in Section 3, the Dirichlet representation 11 naturally produces a boundary term with weight \(w_\partial\) and is therefore the appropriate setting whenever the regularization should discourage features on \(\partial\Omega\), i.e., \(\beta > 0\). The Neumann representation 8 , by contrast, contains no boundary term, and is the appropriate setting when boundary features are admissible, i.e., \(\beta = 0\). We therefore use the Dirichlet-based weights with \(\beta > 0\) for the inverse source problem and the Neumann-based weights with \(\beta = 0\) for the gravimetry problem.

If not explicitly stated otherwise, we employed \(p=2\) in the definition of the weight functions \(w\) and \(\widetilde{w}\), cf. 2 and 4 .

Inverse Source Problem↩︎

The first simulations concern the recovery of spherical shapes of different radius. We observe in Figure 2 that a source with small radius is quite well recovered when located close to the boundary, and that relatively larger sources are also recoverable deeper into the domain. Let us mention that when we ran experiments with deep sources with small radii, the spatial extent of the reconstructed source was too large and the magnitude too small.

These observations can be linked to assumptions ?? and ?? . Admittedly, it is unrealistic that ?? or ?? is perfectly satisfied in a practical setting involving elliptic PDEs, but they can nevertheless shed some light on recoverability: Panels (a) and (b) of Figure 3 show \((K\!\left[\nabla G(\cdot;x)\cdot {\boldsymbol{\nu}}^*(x)\right])(z_R)\) and \((K\!\left[\nabla G(\cdot;x)\cdot {\boldsymbol{\nu}}^*(x) \right])(z_L)\) as functions of the polar angle of \(x\) located at a small circle (a) and a large circle (b). Here, \(z_R\) (red) and \(z_L\) (blue) are boundary points. In neither of the cases, \((K\!\left[\nabla G(\cdot;x)\cdot {\boldsymbol{\nu}}^*(x)\right])(z_R) \geq 0\) for all \(x\) at the circle or \((K\!\left[\nabla G(\cdot;x)\cdot {\boldsymbol{\nu}}^*(x) \right])(z_L) \geq 0\) for all \(x\) at the circle, but for the large source this is closer to being satisfied.

In Figure 4, we investigate how well different shapes can be recovered by solving 1 or 3 . For the true convex shapes in panels (a) and (d), we observe that the directionally weighted method (seen in panels (c) and (f)) is successful in reconstructing the shape and position, whereas the isotropic functional (seen in panels (b) and (e)) struggles to preserve the sharp boundary.

For the non-convex shapes in panels (g) and (j), we observe that both the isotropic and directionally weighted procedures can recover the position and partially the spatial extent of the sources, but they both fail to reconstruct the exact shapes - and the isotropic method still produces smoother reconstructions.

We also investigated how choosing of \(p = 1\), instead of \(p = 2\), in the definition 4 of \(\widetilde{w}\) affected recovery. In Figure 5 we present the recovery of the same sources as in Figure 4. The position and extension of the sources are still quite well preserved, but we might observe a bit more axis-alignment in the inverse solutions.

Gravimetry inspired problem↩︎

The final experiments still concerns the PDE in 6 , but we are now searching for a discontinuous source extending throughout the entire domain \(\Omega\), placing it into the context of a gravimetry type of problem.

In Figure 6, panel (a) shows the true source. When the solution is penalized at the boundary, i.e., when \(\beta > 0\), the reconstruction degrades, as seen in panels (b) and (d): since the true source genuinely reaches the boundary, the boundary penalty suppresses a feature that should be present, and this happens whether or not spatial weighting is used. With \(\beta = 0\), by contrast, the reconstructions improve, as shown in panels (c), (e), and (f). Although, without spatial weighting (panel (c)), the interface is estimated too close to the boundary. Spatial weighting removes this bias: both the Dirichlet-based weights (panel (e)) and the Neumann-based weights (panel (f)) yield good reconstructions.

However, when we restrict the observation domain to the subset \(\Gamma = \{y \in \partial\Omega: y = 1\} \subset \partial\Omega\), we can observe a significant distinction between the reconstructions depending on whether we use Dirichlet or Neumann conditions for generating the Green’s functions, with the Neumann approach clearly being superior, cf. Figure 7.

6 Conclusion↩︎

In this paper we have studied weighted total variation regularization for linear inverse problems whose forward operator \(K \colon L^2(\Omega) \to L^2(\partial\Omega)\) has a significant null space. We have proposed and analyzed the directionally weighted variant 3 4 , in which the worst-case supremum over directions defining \(w\), see 1 2 , is replaced by the actual jump direction of the candidate solution, yielding a sharper but still convex regularizer \(\operatorname{TV}_{\widetilde{w}}\).

In comparison to standard unweighted TV, we observed large improvements in recoverability for both of the weighted regularization techniques - in particular with respect to the center of mass and spatial extension. However, the directional weighting method did seem to produce better results regarding reconstruction of (some of) the shapes, compared with the use of \(w\).

Both weights are constructed directly from the operator \(K\) together with the Green’s function for the Laplace operator on \(\Omega\): depending on whether this Green’s function is taken to satisfy homogeneous Neumann or homogeneous Dirichlet boundary conditions, the source representation produces either a purely interior term, corresponding to the choice \(\beta = 0\) in the variational problem, or a representation with an additional boundary term coming equipped with a natural boundary weight \(w_\partial\), which corresponds to the choice \(\beta > 0\). Choosing the Green’s function according to whether one wants to penalize boundary jumps therefore aligns the regularizer with the (expected) geometry of the unknown source.

When the support of the true source extends to the boundary, we observed in numerical experiments that using Neumann boundary conditions outperformed the use of Dirichlet boundary conditions. Moreover, for all cases where the true support was strictly interior, we observed that the size and position of all kinds of sources were well recovered, but the new methods admittedly recovered circular/oval shapes much better than more involved shapes.

Although it appears that both the total mass and the center of mass of the true sources are well recovered in our numerical simulations, it is currently an open problem whether we can prove that these quantities can be recovered for some typical shapes – or alternatively, that some error estimates can be given.

Also, if more prior information concerning the shape(s) of the potential sources is known, it is currently unclear how this can be incorporated into the regularizer – potentially through a more sophisticated choice of a PDE for generating the involved Green’s functions.

AI Declaration↩︎

The software used in this project was developed with the aid of AI.

References↩︎