dot?

# Vector-valued Spline Method for the Spherical Multiple-shell Electro-
magnetoencephalography Problem

S. Leweke  ,  O. Hauk  and  V. Michel  Geomathematics Group, Department of
Mathematics, University of Siegen, Walter-Flex Straße 3, 57068 Siegen  MRC
Cognition and Brain Sciences Unit, University of Cambridge, 15 Chaucer Road,
Cambridge CB2 7EF, UK  Geomathematics Group, Department of Mathematics,
University of Siegen, Walter-Flex Straße 3, 57068 Siegen
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Human brain activity is based on electrochemical processes, which can only be
measured invasively. For this reason, quantities such as magnetic flux density
(via MEG) or electric potential differences (via EEG) are measured non-
invasively in medicine and research. The reconstruction of the neuronal
current from the measurements is a severely ill-posed problem though the
visualization of the cerebral activity is one of the main research tools in
cognitive neuroscience. Here, using an isotropic multiple-shell model for the
geometry of the human head and a quasi-static approach for modeling the
electro-magnetic processes, we derive a novel vector-valued spline method
based on reproducing kernel Hilbert spaces (RKHS) in order to reconstruct the
neuronal current from the measurements. The presented vector spline method
follows the path of former spline approaches and provides classical minimum
norm properties. In addition, it minimizes the (infinite-dimensional)
Tikhonov-Philips functional which handles the instability of the inverse
problem. This optimization problem reduces to solving a finite-dimensional
system of linear equations without loss of information, due to its particular
construction. It results in a unique solution which takes into account that
only the harmonic and solenoidal component of the neuronal current affects the
measurements. In addition, we prove a convergence result: the solution
achieved by the vector spline method converges to the generator of the data as
the number of measurements increases. The vector splines are applied to the
inversion of several synthetic test cases, where the irregularly distributed
data situation could be handled very well. Combined with several parameter
choice methods, numerical results are shown for synthetic test cases with and
without additional Gaussian white noise. Former approaches based on scalar
splines are outperformed by the novel vector splines results with respect to
the normalized root mean square error. Finally, results for real data are
demonstrated. They can be computed quickly and are reasonable with respect to
physiological expectations.

###### 2020 Mathematics Subject Classification:

41A15, 42C10, 45B05, 46C07, 46N40, 47A52, 65D07, 65R30, 65R32

Keywords. electroencephalography, ill-posed problems, integral equation,
inverse problems, magnetoencephalography, regularization methods, reproducing
kernel Hilbert spaces, scalar spherical splines, vector spherical splines

##  1\. Introduction

Neurons in the cerebrum use electrochemical processes for signaling. These
electrical activities are often referred to as cerebral or brain activity,
respectively. First of all, the signal processing and interplay of neuronal
cells in the brain is responsible for human behavior and cognition.
Accordingly, defects in the signal processing cause diseases, such as epilepsy
or schizophrenia. Second, localization of the brain activity can help to
answer open questions in the field of cartography of the brain structures and
to further understand the human brain functionality. Hence, an accurate
visualization of the cerebral activity is one of the main tools in research
and diagnosis, [33].

Since non-invasive measurements outside the head are often the only available
option, the physical fundamentals connecting the brain activity with
measurable quantities need to be understood and modeled adequately. The
mathematical modeling of these fundamentals results in an electric current
inside the cerebrum which can be split into the neuronal current $J$ and its
induced Ohmic current, where only the neuronal current is of interest [32].
The cerebral current induces a macroscopic electric potential $u_{s}$ and a
magnetic field $B$, which are transmitted through the conductive brain tissues
and can be measured outside the head non-invasively provided that the number
of simultaneously active neurons is large enough.

In order to formulate the relation between the neuronal current and the
measured quantities, we follow the path of optimize-then-discretize approaches
in [35, 34, 13, 19, 20, 41]. Via this ansatz, the problem is solved
analytically as far as possible. This has the advantage that structures like
non-visible parts of the current are conserved up to the implementation.
Otherwise, model assumptions are required to enable the analysis. For our
approach, it is sufficient to model the neuronal current as continuously
distributed current, the head via the common multiple-shell model, [30, 41],
and the physics via quasi-static Maxwell’s equations, [30, 45, 46].

Within this framework, the problem of directly reconstructing the vector-
valued neuronal current from magnetic flux density $\nu\cdot B$ measurements
by the magnetoencephalograph (MEG), [11], and electric potential differences
values on the scalp obtained via the electroencephalograph (EEG), [12], was
addressed in [35, 34]. Therein, singular value decompositions (SVD) of the
corresponding operators are derived based on a novel set of vector-valued
orthonormal basis functions. First of all, the SVD paves the path for
determining the parts of the neuronal current which can be reconstructed from
the measurements. This is closely related to the specification of the operator
null space. Moreover, asymptotic behavior of the singular values characterizes
the ill-posedness of the inverse problems. Related problems of reconstructing
scalar parts of the cerebral current from these measurements have been
discussed before in [13, 19, 20, 41].

Based on the insights achieved by SVDs, it is revealed that for the
reconstruction of the neuronal current $J$ from the given measurements several
challenges need to be tackled:

Non-uniqueness:

If a solution of the inverse problem exists, it is not unique, since only the
harmonic part of the solenoidal directions of the neuronal current is not
silent for the MEG and EEG, [34].

Instability:

Due to the distance between the head and the MEG device and the distance
between the cerebrum and the outer surface of the scalp, terms emerge in the
formulae which are associated to a continuation of the magnetic/electric field
down to a sphere with lower radius. This is a well-known mathematical problem,
which also occurs in the geosciences for the magnetic and the gravitational
field. It leads to exponentially diverging singular values for the inverse
operator, which implies a severe ill-posedness of the inverse problem. As a
consequence, small changes in the data (e.g. by means of noise) have a huge
impact on the reconstruction.

Noise:

The achieved data is noisy, due to physiological noise (e.g. head, eye, jaw,
or neck movements, magnetic field generated by the heart) and non-
physiological noise (e.g. technical noise level of the devices, particles
attached or implanted to the patient’s body, other outer fields and currents),
[30, 34].

Data Distribution:

There are only few sensor positions available (commonly a few dozens to a few
hundreds in the case of only MEG or EEG measurements), which are irregularly
distributed and have a major gap in the area of the face and the lower half of
the head.

The aim of this paper is to reconstruct the vector-valued neuronal current
directly by developing a robust and stable numerical optimize-then-discretize
method for solving these functional inverse problems addressing these
challenges. Within this paper, we will mainly answer the following questions

  * •

Does our method produce a reasonable and correct reconstruction?

  * •

Is the method stable with respect to increasing noise level?

  * •

Is there an advantage of reconstructing the vector-valued current directly
instead of reconstructing scalar-valued components as an intermediate step?

First of all, regularization techniques are required to handle the ill-
posedness: the best-approximate solution obtained from noisy data is not
suitable for the reconstruction, since the generalized Moore-Penrose inverse
is unbounded if it exists, [18]. Popular regularizations include generalized
Tikhonov-Philips functionals. In the context of approximations of the neuronal
current based on continuously distributed MEG and EEG models, a regularization
method based on global orthonormal basis functions (related to spherical
harmonics) and scalar spherical splines have been used before, [20]. In
addition, for the EEG model a hybrid analytical-numerical algorithm using
OpenMEG exists, [31]. Besides, discrete models (e.g. MNE, LORETA, FOCUSS)
often yield in a finite-dimensional Tikhonov-regularized normal equation,
which is solved by a variety of optimization algorithms, [28, 32]. Though such
methods have also become established, their usage has to be seen critically
from the mathematical point of view, since the discretization approach does
not consider the known results on the large null spaces of the inverse
problems. Hence, it is not possible to distinguish uniquely determined
components of the solution from remaining degrees of freedom including
possible artifacts.

Within this paper, we extend the idea of the regularization method based on
global orthonormal basis functions and construct vector-valued as well as
scalar-valued splines based on the SVD of the integral operators. In order to
obtain a unique solution, we assume additionally that the neuronal current
satisfies an $\mathrm{L}^{2}$-minimum-norm condition. This results (in the
unregularized case) in the best-approximate solution of the inverse problems,
[34, Thm. 13.6], which is entirely contained in the orthogonal complement of
the operator null space.

The idea of using a reproducing kernel to regularize an inverse problem is not
new. For instance, [43] presented a regularization theory for linear inverse
problems in RKHSs. Regarding examples of approaches tailored for specific
domains, (functional) spherical splines based on reproducing kernels (which
will shortly be referred to _spherical splines_) are a widely known tool for
interpolation and approximation, which is frequently used for solving
(inverse) problems in the geosciences and beyond. This method goes back to
[23] and [50] and provides us with a best-approximate solution, [23, 26, 24,
50], in the unregularized case. They are constructed via spatially localized
reproducing kernels, hence local changes in the data have only mainly local
effects in the reconstruction. In addition, spherical splines have the
smallest (RKHS) norm among all functions fitting the data. In addition, the
formulation via splines allows, without loss of information, to transfer the
infinite-dimensional optimization problem of minimizing the Tikhonov-Philips
functional to a finite-dimensional system of linear equations. Due to these
nice properties, scalar spherical splines have been used in several
applications on the unit sphere, [39, 38, 25, 48], and on the ball, [2, 8, 9].
Also in the MEG and EEG setting they have been used before in [20].

In this paper, we will construct mainly two types of spherical splines over
the three-dimensional ball. For the first type, further decompositions of the
neuronal current are required in order to achieve a reduction to a scalar-
valued problem. Then, scalar-valued splines connected to tailor-made
orthonormal basis functions are constructed to solve these problems. This
follows the path of [20] but it takes into account the novel knowledge about
those parts of the neuronal current which are silent to MEG and EEG.
Afterwards, the scalar solutions have to be transferred back to the neuronal
current which we are interested in. The second and novel type of spherical
splines on the ball is vector-valued and can be used to solve the vector-
valued MEG and EEG problem directly. This novel approach conserves the
excellent properties of former (scalar) spline methods like the minimum
properties.

From a theoretical point of view, a direct vector approach for the entire
neuronal current has several advantages over the scalar approaches. First of
all, one obtains a direct relation between the desired quantity and the
measured one which reduces the effort of transforming the scalar quantities
back to the current. This is an advantage, since in several cases, there
exists no one-to-one relation between the scalar and the vector quantities.
This takes effects especially when adding uniqueness constraints. Besides,
additional smoothness or boundary conditions or complementary gauges, which
may not have a physical or medical meaning, are required for the current in
order to derive or transfer the scalar solutions. Finally, from the vector
approach and its characterization of the null space it became explicit that
only the combined inversion of MEG and EEG data can depict a comprehensive
image of the current. The vector approach yields the easiest possibility to
combine these to inversions compared to the scalar approaches, see [34, 35].

Within this paper, we will show that also from a numerical point of view, the
vector approach has several advantages compared to the scalar ones which
manifests in lower approximation errors and a higher stability with respect to
the influence of noise. This is done by several synthetic numerical tests.
Afterwards the vector spline method is applied to real data, where the
corresponding approximation of the neuronal current answers all physiological
expectations.

The rest of the paper is structured as follows: In Sec. 2 short introductions
into (vector-valued) orthonormal basis functions with generalized Fourier
analysis on the ball, the multiple-shell model for the MEG and EEG problem,
and the scalar reproducing kernel based spline method for inverse problems are
given. In Sec. 3 vector-valued functional RKHS and the corresponding splines
on the ball are constructed. Afterwards, properties of the spline such as
minimum properties, approximation property, and convergence results are
proven. In Sec. 4 the foundations for the implementation of the particular MEG
and EEG scalar and vector spline methods are summarized. The particular
numerical framework and the numerical results for the synthetic test cases as
well as the inversion of the real data are stated in Sec. 5. Therein, the used
methods are also checked against each other. Eventually, conclusions are
presented in Sec. 6.

##  2\. Preliminaries

###  2.1. Construction of Suitable Orthonormal Basis Functions

As stated before, an adequate set of orthonormal basis functions is required
for the SVD of the integral operators. In addition, due to Mercers
representation theorem, [37], it is well-known that reproducing kernels are
closely connected to orthonormal basis functions. Eventually, these particular
basis functions build the foundation of the presented spline methods. Before
we have a closer look at these basis functions, we give a short introduction
to the used notations.

The Euclidean inner product is denoted by $\cdot$, the vector product by
$\times$, and the tensor product by $\otimes$. In addition, we will use
throughout the paper the abbreviations $x=r\xi$ and $y=s\eta$ with unit
vectors $\xi$, $\eta\in S\coloneqq S_{1}\subset\mathbb{R}^{3}$ and radii
$r=|x|$, $s=|y|$, where $S_{1}$ is the unit sphere in $\mathbb{R}^{3}$. This
decomposition is unique for all $x\in\mathbb{R}^{3}\setminus\\{0\\}$.
Furthermore, $B_{R}\subset\mathbb{R}^{3}$ is the ball with radius $R$ and
center $0$. In our notation, $\nabla^{\ast}_{\xi}$ denotes the part of the
gradient $\nabla$ containing the tangential derivatives divided by $r$, [27,
Eq. (2.136)]. $\nabla^{\ast}_{\xi}$ is often called the surface gradient. The
differential operator $L^{\ast}$ (independent of the radius $r$) is defined by

$L^{\ast}_{\xi}\coloneqq x\times\nabla_{x}=\xi\times\nabla^{\ast}_{\xi}$ (2.1)

and called the surface curl operator. The Beltrami operator is given by
$\Delta^{\ast}_{\xi}=\nabla^{\ast}_{\xi}\cdot\nabla^{\ast}_{\xi}$, [27, Eq.
(2.140)], and is the part of the Laplacian independent of the radius. Note
that variables as indices of operators indicate the dependence to which the
operator is implied, since this is otherwise not always unique (e.g. (2.4)
below).

The construction of our basis function is in all cases based on spherical
harmonics. Recall that a function of $\mathrm{Harm}_{n}(S)$ (i.e. the space of
all homogeneous, harmonic polynomials of degree $n$ restricted to the unit
sphere $S$, [27, Def. 3.22]) is called a spherical harmonic of degree
$n\in\mathbb{N}_{0}\coloneqq\mathbb{N}\cup\mathchoice{\left\\{0\right\\}}{\left.\\{0\\}\right.}{\left.\\{0\\}\right.}{\left.\\{0\\}\right.}$.
With
$\mathchoice{\left\\{Y_{n,j}\right\\}_{j=1,\dots,2n+1}}{\left.\\{Y_{n,j}\\}_{j=1,\dots,2n+1}\right.}{\left.\\{Y_{n,j}\\}_{j=1,\dots,2n+1}\right.}{\left.\\{Y_{n,j}\\}_{j=1,\dots,2n+1}\right.}$
we denote an $\mathrm{L}^{2}(S)$-orthonormal set in $\mathrm{Harm}_{n}(S)$,
[27, Rem. 3.25]. For more details on scalar spherical harmonics and their
properties, see for instance [27]. For an introduction to (vector-valued)
Lebesgue spaces, see [5].

By means of the scalar spherical harmonics, we can define a complete
$\mathrm{L}^{2}(S,\mathbb{R}^{3})$-orthonormal system of vector-valued
spherical harmonics, [27, Thm. 5.56], which goes back to Edmonds, [16], and is
defined for example in [27, Eq. (5.309)-(5.311)] by

$\tilde{{y}}_{n,j}^{(i)}(\xi)\coloneqq\left(\tilde{\mu}_{n}^{(i)}\right)^{-1/2}\tilde{o}^{(i)}_{n,\xi}Y_{n,j}(\xi),\qquad\xi\in
S,$ (2.2)

where

$\tilde{\mu}_{n}^{(i)}\coloneqq\begin{cases}(n+1)(2n+1),&\text{for }i=1,\\\
n(2n+1),&\text{for }i=2,\\\ n(n+1),&\text{for
}i=3,\end{cases}\qquad\tilde{o}^{(i)}_{n,\xi}\coloneqq\begin{cases}(n+1)\xi-\nabla^{\ast}_{\xi},&\text{for
}i=1,\\\ n\xi+\nabla^{\ast}_{\xi},&\text{for }i=2,\\\
L^{\ast}_{\xi},&\text{for }i=3,\end{cases}$ (2.3)

for all $i\in\\{1,2,3\\}$, $n\in\mathbb{N}_{0_{i}}$, and $j=1,\dots,2n+1$. We
use $\mathbb{N}_{0_{i}}$ as an abbreviation for $\mathbb{N}_{0}$ in the case
of $i=1$ and for $\mathbb{N}$ in the case of $i\in\\{2,3\\}$. Note that these
vector-valued spherical harmonics are homogeneous harmonic polynomials. This
and more information on vector spherical harmonics can be found in [27]. In
analogy, we can define, for all $i\in\\{1,2,3\\}$ and
$n\in\mathbb{N}_{0_{i}}$, Edmonds-vector-Legendre polynomials by means of the
scalar Legendre polynomials $P_{n}$, [27, Eq. (3.165)], of degree
$n\in\mathbb{N}_{0}$ and type $i\in\\{1,2,3\\}$ as in [27, Lem. 5.63],

$\tilde{{p}}^{(i)}_{n}({\xi},{\eta})\coloneqq\left(\tilde{\mu}_{n}^{(i)}\right)^{-1/2}\tilde{o}^{(i)}_{n,\xi}P_{n}({\xi}\cdot{\eta}),\qquad\xi,\,\eta\in
S.$ (2.4)

Besides this, the Legendre polynomials and their vectorial counterpart enable
addition theorems, [27, Thm. 3.26, Thm. 5.64]. These imply for all types and
degrees the representations

$\sum_{j=1}^{2n+1}Y_{n,j}({\xi})Y_{n,j}({\eta})=\frac{2n+1}{4\pi}P_{n}({\xi}\cdot{\eta}),\qquad\sum_{j=1}^{2n+1}\tilde{{y}}^{(i)}_{n,j}({\xi})Y_{n,j}({\eta})=\frac{2n+1}{4\pi}\tilde{{p}}^{(i)}_{n}({\xi},{\eta}),\qquad\xi,\,\eta\in
S.$ (2.5)

Furthermore, we also obtain for (pointwise) Euclidean norms

$\sum_{j=1}^{2n+1}\mathchoice{\left\lvert\tilde{{y}}^{(i)}_{n,j}({\xi})\right\rvert^{2}}{\lvert\tilde{{y}}^{(i)}_{n,j}({\xi})\rvert^{2}}{\lvert\tilde{{y}}^{(i)}_{n,j}({\xi})\rvert^{2}}{\lvert\tilde{{y}}^{(i)}_{n,j}({\xi})\rvert^{2}}=\frac{2n+1}{4\pi},\qquad\xi\in
S.$ (2.6)

Based on a separation ansatz, we combine the vector spherical harmonics with
orthogonal Jacobi functions $P_{m}^{(\alpha,\beta)}$ for the radial part in
order to achieve an orthonormal set of functions over the ball. For further
details on Jacobi polynomials, see, for instance, [49].

######  Theorem 2.1.

The set of vector-valued functions $\tilde{{g}}^{(i)}_{m,n,j}(R;\cdot)$ for
$i\in\\{1,2,3\\}$, $m\in\mathbb{N}_{0}$, $n\in\mathbb{N}_{0_{i}}$, and
$j=1,\dots,2n+1$ with the parameter

$t_{n}^{(i)}\coloneqq\begin{cases}n&i=1,\,3,\\\ n-1&i=2,\\\ \end{cases}$

is defined via

$\tilde{{g}}^{(i)}_{m,n,j}(R;{x})\coloneqq\sqrt{\frac{4m+2t^{(i)}_{n}+3}{R^{3}}}\left(\frac{r}{R}\right)^{t^{(i)}_{n}}P_{m}^{\left(0,t^{(i)}_{n}+1/2\right)}\left(2\frac{r^{2}}{R^{2}}-1\right)\tilde{{y}}_{n,j}^{(i)}(\xi),\qquad
x\in B_{R}.$

It is a complete orthonormal system in $\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})$.
Eventually, each $f\in\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})$ has with the
abbreviation $f^{\wedge}(i,m,n,j)\coloneqq\left\langle
f,\tilde{{g}}^{(i)}_{m,n,j}(R;\cdot)\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}$
the representation

$f=\sum_{i=1}^{3}\sum_{m=0}^{\infty}\sum_{n=0_{i}}^{\infty}\sum_{j=1}^{2n+1}f^{\wedge}(i,m,n,j)\tilde{{g}}^{(i)}_{m,n,j}(R;\cdot),$
(2.7)

which converges unconditionally and strongly in the
$\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})$-sense.

Recall that other choices for the parameter sequence $(t_{n}^{(i)})_{n}$ are
possible, as long as the condition
$\inf_{n\in\mathbb{N}_{0_{i}}}t_{n}^{(i)}\geq-\frac{3}{2}$ is satisfied, see
[34], also regarding a proof of Thm. 2.1. However, the stated sequence is the
natural choice corresponding to our application, [35, Rem. 5.2].

######  Definition 2.2.

Let $f\in\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})$ be a given function. For almost
all $r\in[0,R]$, $i\in\\{1,2,3\\}$, $n\in\mathbb{N}_{0_{i}}$, and
$j=1,\dots,2n+1$ we define

$f_{n,j}^{(i)}(r)\coloneqq\int_{S}f(r\xi)\cdot\tilde{y}_{n,j}^{(i)}(\xi)\,\mathrm{d}\omega(\xi).$
(2.8)

An analogue can be obtained for scalar-valued functions
$F\in\mathrm{L}^{2}(B_{R})$ with an appropriate sequence
$(t_{n})_{n\in\mathbb{N}_{0}}$ by

$F_{n,j}(r)\coloneqq\int_{S}F(x)Y_{n,j}(\xi)\,\mathrm{d}\omega(\xi).$ (2.9)

###  2.2. The Multiple-shell Model for MEG and EEG

Having the described scalar- and vector-valued orthonormal basis function at
hand, we are able to recapitulate the SVD of the integral operators and
related relations between the measured quantities (i.e., magnetic flux density
and electric potential differences) and the neuronal current $J$. This SVD is
a central point for the construction of the (vector) functional splines for
inverse problems.

Before, we state these results, we introduce the common _multiple-shell model_
used for modeling the head, [42]. More precisely, in our setting, we assume
that

  * •

the cerebrum is a closed ball with radius $\varrho_{0}$, that is
$B_{\varrho_{0}}$,

  * •

surrounding the cerebrum, there are $L\geq 2$ spherical shells
$S_{[\varrho_{l},\varrho_{l+1}]}\coloneqq\overline{B_{\varrho_{l+1}}\setminus
B_{\varrho_{l}}}$ for $l=1,\dots,L-1$ modeling the various head tissues,

  * •

each tissue (i.e. each shell $S_{[\varrho_{l},\varrho_{l+1}]}$) has a constant
conductivity $\sigma_{l}$ for all $l=0,\dots,L-1$ and outside the head the
conductivity is vanishing, $\sigma_{L}=0$,

  * •

the permeability is constant everywhere and equals the permeability of the
vacuum $\mu_{0}$,

  * •

the relation between the neuronal current and the induced quantities can be
modeled by means of quasi-static Maxwell’s equations, [45], and

  * •

the continuously distributed neuronal current
$J\in\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})$ is non-vanishing only
inside the cerebrum.

Note that we do not assume further smoothness or boundary conditions for the
neuronal current.

In the MEG case, the functionals
$\mathcal{A}_{\mathrm{M}}^{k}\colon\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})\to\mathbb{R}$
map the neuronal current onto the magnetic flux density evaluated at the
sensor positions $y_{k}\in\overline{\mathbb{R}^{3}\setminus B_{\varrho_{L}}}$,
$k=1,\dots,\ell_{M}$, outside the head. The flux density is the part of the
magnetic field pointing towards the normal direction $\nu$ of the sensor
surface. In the EEG case, the functionals
$\mathcal{A}_{\mathrm{E}}^{k}\colon\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})\to\mathbb{R}$,
$k=1,\dots,\ell_{\mathrm{E}}$, map the current onto the electric potential
difference measured at several positions on the scalp.

Summarizing the results, we obtain the following equations for functionals,
where they are stated in previous publications, [34, Eq. (6.5)] in the MEG
case and in [34, Proof of Thm. 6.3] in the EEG case, respectively:

$\displaystyle\mathcal{A}_{\mathrm{M}}^{k}J$ $\displaystyle=\nu(y_{k})\cdot
B(y_{k})$
$\displaystyle=-\mu_{0}\sum_{n=1}^{\infty}\sum_{j=1}^{2n+1}\sqrt{\frac{n\varrho_{0}}{(2n+1)(2n+3)}}\mathchoice{{\left\langle
J,\tilde{g}^{(3)}_{0,n,j}(\varrho_{0};\cdot)\right\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}^{(3)}_{0,n,j}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}^{(3)}_{0,n,j}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}^{(3)}_{0,n,j}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}\left(\frac{\varrho_{0}}{s_{k}}\right)^{n+1}\frac{1}{s_{k}}\nu(y_{k})\cdot\tilde{y}_{n,j}^{(1)}(\eta_{k}),$
(2.10) $\displaystyle\mathcal{A}_{\mathrm{E}}^{k}J$
$\displaystyle=u_{s}(y_{k})$
$\displaystyle=\sum_{n=1}^{\infty}\sum_{j=1}^{2n+1}\sqrt{\frac{n}{\varrho_{0}}}\mathchoice{{\left\langle
J,\tilde{g}_{0,n,j}^{(2)}(\varrho_{0};\cdot)\right\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}_{0,n,j}^{(2)}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}_{0,n,j}^{(2)}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}_{0,n,j}^{(2)}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}\left(\frac{(n+1)}{n}\left(\frac{s_{k}}{\varrho_{L}}\right)^{2n+1}+1\right)\left(\frac{\varrho_{0}}{s_{k}}\right)^{n+1}\beta^{(L)}_{n}Y_{n,j}(\eta_{k}).$
(2.11)

Recall that $y_{k}=s_{k}\eta_{k}\in\mathbb{R}^{3}$ for all
$k=1,\dots,\ell_{\mathrm{M}}$ or $k=1,\dots,\ell_{\mathrm{E}}$, respectively,
denotes the sensor positions of the measurement devices. Note that the
coefficients $(\beta^{(L)}_{n})_{n}$ depend on the particular model geometry
and the tissues conductivities. A recursive formula for determining these
coefficients and an analysis of the corresponding asymptotic behavior are
given in [34, 35]. Therein, the absolute and uniform convergence of the series
are additionally proven.

Before introducing the corresponding scalar-valued problems, we shortly sum up
the information which can be achieved by the SVD-based representation of the
functionals. In the MEG case, (2.10) reveals that only the direction of the
neuronal current which corresponds to the toroidal part (i.e $i=3$) is not
silent for the MEG device. In complementary, in the EEG case, (2.11) shows
that only the orthonormal basis functions of type $i=2$ affect the
measurements. Having the construction of Edmonds vector spherical harmonics in
mind and the fact that only type $i=2$ and $i=3$ are divergence-free vector
fields, one can deduce that only the solenoidal direction of the neuronal
current is related to the two measurements. In addition, not the entire
solenoidal part of the current can be reconstructed, since only the degree
$m=0$ of the generalized Fourier expansions contribute to the data. This
coincides with the harmonic parts of the non-silent directions. Hence, the
native non-uniqueness condition is to require a harmonic and solenoidal
neuronal current. This also coincides with the minimum-norm condition, which
will serve as our uniqueness constraint for the neuronal current. In addition,
the severe ill-posedness of the two inverse problems can be seen in the
exponentially fast decreasing of the singular values to zero, which mainly
goes back to the sequence $((\varrho_{0}/s_{k})^{n+1})_{n}$. A more detailed
discussion of the non-unique solution, additional uniqueness constraints, and
the singular values can be found in [35, 34].

Besides a decomposition of the neuronal current by means of this vector-valued
orthonormal basis, several other decompositions exist with the aim of only
considering relevant scalar-valued parts of the neuronal current. For example,
the Helmholtz decomposition with the Coulomb gauge of the neuronal current can
be used for the MEG as well as the EEG problem, [14, 20, 21, 19]. Note that
for the MEG problem, in [20], the Helmholtz decomposition is combined with a
layer density constraint to achieve uniqueness of the solution. Unfortunately,
this uniqueness assumption contradicts our minimum norm assumption, [35, Sec.
20.2]. Thus, we need to adapt the approach in order to fit into our setting.

In general, the Helmholtz decomposition for a sufficiently smooth vector-
valued function is given by

${J}={\nabla}\Psi+{\nabla}\times{a}$ (2.12)

with the scalar potential $\Psi\in\mathrm{C}^{2}(B_{\varrho_{0}})$ and the
vector potential $a\in\mathrm{C}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})$. Note
that this decomposition is not unique without an additional gauge, therefore,
we use the Coulomb gauge (i.e. $\nabla\cdot a=0$). Due to its smoothness, the
vector potential $a$ is decomposable further by means of the spherical
Helmholtz decomposition, [27, Eq. (5.58)],

$a(x)=\xi A^{(1)}(x)+\nabla^{\ast}_{\xi}A^{(2)}(x)+L^{\ast}_{\xi}A^{(3)}(x).$

Inserting the decompositions for the neuronal current into the functionals
(2.10), (2.11) and using several orthogonality properties of the orthonormal
basis function, [35, 34], we obtain new functionals derived from the previous
ones. In the MEG case, we achieve a functional
$\mathcal{A}^{k}_{\mathrm{m}}\colon\mathrm{C}^{2}(B_{\varrho_{0}})\to\mathbb{R}$,
$k=1,\dots,\ell_{\mathrm{M}}$, mapping the _scalar-valued part_ $A^{(1)}$ of
the neuronal current onto the magnetic flux density evaluated at the sensor
positions [34, Lem. 7.3], that is

$\begin{multlined}\mathcal{A}^{k}_{\mathrm{m}}A^{(1)}\coloneqq\mu_{0}\nu(y_{k})\cdot\left(\int_{B_{\varrho_{0}}}(\Delta_{x}(|x|A^{(1)}(x)))\nabla_{y}K_{\mathrm{m}}(x,y)\,\mathrm{d}x\middle)\right|_{y=y_{k}}\\\
=-\mu_{0}\sum_{n=1}^{\infty}\sum_{j=1}^{2n+1}\sqrt{\frac{1}{(2n+1)(n+1)}}\left(\frac{\mathrm{d}A_{n,j}^{(1)}}{\mathrm{d}r}(\varrho_{0})\varrho_{0}-(n-1)A_{n,j}^{(1)}(\varrho_{0})\right)\left(\frac{\varrho_{0}}{s_{k}}\right)^{n+2}\nu(y_{k})\cdot\tilde{y}_{n,j}^{(1)}(\eta_{k}).\end{multlined}\mathcal{A}^{k}_{\mathrm{m}}A^{(1)}\coloneqq\mu_{0}\nu(y_{k})\cdot\left(\int_{B_{\varrho_{0}}}(\Delta_{x}(|x|A^{(1)}(x)))\nabla_{y}K_{\mathrm{m}}(x,y)\,\mathrm{d}x\middle)\right|_{y=y_{k}}\\\
=-\mu_{0}\sum_{n=1}^{\infty}\sum_{j=1}^{2n+1}\sqrt{\frac{1}{(2n+1)(n+1)}}\left(\frac{\mathrm{d}A_{n,j}^{(1)}}{\mathrm{d}r}(\varrho_{0})\varrho_{0}-(n-1)A_{n,j}^{(1)}(\varrho_{0})\right)\left(\frac{\varrho_{0}}{s_{k}}\right)^{n+2}\nu(y_{k})\cdot\tilde{y}_{n,j}^{(1)}(\eta_{k}).$
(2.13)

The precise representation of the integral kernel $K_{\mathrm{m}}$ can be
found in (5.5).

Apart from that, inserting the Helmholtz decomposition into the functionals
for the EEG problem, two scalar parts of the neuronal current affect the
measured quantity. On the one hand, we have the scalar potential $\Psi$ of the
Helmholtz decomposition and on the other hand, we have the scalar part
$A^{(3)}$ of the vector potential. The functionals
$\mathcal{A}^{k}_{\mathrm{e}}\colon\mathrm{C}^{2}(B_{\varrho_{0}})\times\mathrm{C}^{2}(B_{\varrho_{0}})\to\mathbb{R}$,
$k=1,\dots,\ell_{\mathrm{e}}$ are given by [34, Sec. 7.6]

$\mathcal{A}^{k}_{\mathrm{e}}\left(\Psi,A^{(3)}\right)=\sum_{n=1}^{\infty}\sum_{j=1}^{2n+1}\left(n\Psi_{n,j}(\varrho_{0})-n(n+1)A^{(3)}_{n,j}(\varrho_{0})\right)\left(\frac{n+1}{n}\left(\frac{s_{k}}{\varrho_{L}}\right)^{2n+1}+1\right)\left(\frac{\varrho_{L}}{s_{k}}\right)^{n+1}\beta_{n}^{(L)}Y_{n,j}(\eta_{k}).$

Note that no gauge is required to derive this representation. In the case of
the Coulomb gauge, only a relation between the functions $A^{(1)}$ and
$A^{(2)}$ can be achieved [35, Thm. 15.14]. In the case of the Poincaré gauge
(i.e. $x\cdot a=0$), only conditions for the function $A^{(1)}$ are gained. In
order to get rid of the function $A^{(3)}$ in the upper relation, an
additional boundary condition for the neuronal current, $(x\cdot
J)|_{S_{\varrho_{0}}}=0$ is often used, to connect the two scalar parts:
$0=\frac{\mathrm{d}}{\mathrm{d}r}\Psi(x)+r^{-1}\Delta^{\ast}_{\xi}A^{(3)}(x)$,
[14, 19, 21]. However, it can be proven that every neuronal current
$J\in\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})$ satisfying this boundary
condition in addition to the minimum norm condition, must be equal to the zero
function, [34, Thm. 15.25]. Hence, the Helmholtz decomposition is not suitable
for combining the MEG and EEG inversion under the minimum-norm condition.

###  2.3. Introduction to Reproducing Kernel Based Spline Methods for Inverse
Problems

The spline method presented in this paper is based on the construction of
reproducing kernels and their generated Hilbert spaces. An overview of
properties of the scalar as well as the vector-valued spherical interpolating
splines based on RKHS on the sphere is summarized in [26].

In the work of [2], the scalar spherical splines are extended to arbitrary
RKHS over compact domains and combined with functionals for the application to
(ill-posed) inverse problems, see [2, 38]. These results are also summarized
for scalar Hilbert spaces over the ball in the particular context of the
inverse MEG and EEG problem in [20].

Before we extend this approach, we briefly summarize the properties and
statements relevant for our application in the following. More information can
be found in the above references and the references therein. Now, let the
system of functions

$H_{n,j}(x)\coloneqq G_{n}(r)Y_{n,j}(\xi),\qquad
n\in\mathbb{N}_{0},\,j=1,\dots,2n+1$

be a linearly independent, orthonormal system in $\mathrm{L}^{2}(B_{R})$. This
is suitable for isotropic applications like the inverse MEG and EEG problem.

Let $\mathscr{H}\subset\mathrm{L}^{2}(B_{R})$ be an RKHS over the ball $B_{R}$
generated by the sequence $(\kappa_{n})_{n}$ and equipped with the inner
product

$\mathchoice{{\left\langle F,G\right\rangle_{\mathscr{H}}}}{{\langle
F,G\rangle_{\mathscr{H}}}}{{\langle F,G\rangle_{\mathscr{H}}}}{{\langle
F,G\rangle_{\mathscr{H}}}}\coloneqq\sum_{n\in\mathbb{N}_{0}}\sum_{j=1}^{2n+1}\kappa_{n}^{2}\mathchoice{{\left\langle
F,H_{n,j}\right\rangle_{\mathrm{L}^{2}(B_{R})}}}{{\langle
F,H_{n,j}\rangle_{\mathrm{L}^{2}(B_{R})}}}{{\langle
F,H_{n,j}\rangle_{\mathrm{L}^{2}(B_{R})}}}{{\langle
F,H_{n,j}\rangle_{\mathrm{L}^{2}(B_{R})}}}\mathchoice{{\left\langle
G,H_{n,j}\right\rangle_{\mathrm{L}^{2}(B_{R})}}}{{\langle
G,H_{n,j}\rangle_{\mathrm{L}^{2}(B_{R})}}}{{\langle
G,H_{n,j}\rangle_{\mathrm{L}^{2}(B_{R})}}}{{\langle
G,H_{n,j}\rangle_{\mathrm{L}^{2}(B_{R})}}}.$

This construction goes back to [24]. The reproducing kernel is uniquely given
via Mercers representation theorem [37] by

$K(x,z)=\sum_{\begin{subarray}{c}n\in\mathbb{N}_{0}\\\ \kappa_{n}\neq
0\end{subarray}}\sum_{j=1}^{2n+1}\kappa_{n}^{-2}H_{n,j}(x)H_{n,j}(z),\qquad
x,\,z\in\mathrm{B}_{R}.$

In order to make sense of the latter series expressions, the sequence
$(\kappa_{n})_{n}$ needs to satisfy a summability condition, that is

$\sum_{\begin{subarray}{c}n\in\mathbb{N}_{0}\\\ \kappa_{n}\neq
0\end{subarray}}\kappa_{n}^{-2}\sup_{x\in
B_{R}}\mathchoice{\left\lvert\sum_{j=1}^{2n+1}H_{n,j}(x)\right\rvert^{2}}{\lvert\sum_{j=1}^{2n+1}H_{n,j}(x)\rvert^{2}}{\lvert\sum_{j=1}^{2n+1}H_{n,j}(x)\rvert^{2}}{\lvert\sum_{j=1}^{2n+1}H_{n,j}(x)\rvert^{2}}=\sum_{\begin{subarray}{c}n\in\mathbb{N}_{0}\\\
\kappa_{n}\neq
0\end{subarray}}\frac{2n+1}{4\pi}\kappa_{n}^{-2}\sup_{r\in[0,R]}\mathchoice{\left\lvert
G_{n}(r)\right\rvert^{2}}{\lvert G_{n}(r)\rvert^{2}}{\lvert
G_{n}(r)\rvert^{2}}{\lvert G_{n}(r)\rvert^{2}}<\infty.$

According to [20, Eq. (41)], the corresponding _(scalar) spline function $S$_
for the interpolation problem

$g_{k}=\mathcal{A}^{k}F,\qquad k=1,\dots,\ell,$

where each linear functional $\mathcal{A}^{k}$ for $k=1,\dots,\ell$ maps from
$\mathscr{H}$ continuously to $\mathbb{R}$, is of the form

$S\coloneqq\sum_{k=1}^{\ell}\alpha_{k}\mathcal{A}^{k}_{z}K(\cdot,z)=\sum_{k=1}^{\ell}\alpha_{k}\sum_{\begin{subarray}{c}n\in\mathbb{N}_{0}\\\
\kappa_{n}\neq
0\end{subarray}}\sum_{j=1}^{2n+1}\kappa_{n}^{-2}\left(\mathcal{A}^{k}H_{n,j}\right)H_{n,j}$
(2.14)

with arbitrary but real coefficients $\alpha=(\alpha_{k})_{k}$. It is well-
known that a scalar spherical spline function has the following properties:

  1. (1)

For given data $g\in\mathbb{R}^{\ell}$ and continuous linear functionals
$\mathcal{A}^{k}$ for $k=1,\dots,\ell$, the interpolation problem is uniquely
solvable if and only if the functionals $\mathcal{A}^{k}$ are linearly
independent, see [2, Thm. 5.10].

  2. (2)

Among all solutions $F\in\mathscr{H}$ interpolating the data, that is
$\mathcal{A}^{k}F=g_{k}$ for $k=1,\dots,\ell$, the spline function is the only
solution with minimal $\mathscr{H}$-norm, see [2, Thm. 10.14].

  3. (3)

The spline function satisfies a best-approximation property, see [2, Thm.
10.16].

In addition, in the regularized case the following theorem holds true.

######  Theorem 2.3 (Spline Approximation, [38, Thm. 10.16]).

Let $g\in\mathbb{R}^{\ell}$ and a regularization parameter $\lambda>0$ be
given. If the vector
$\alpha=(\alpha_{k})_{k=1,\dots,\ell}\in\mathbb{R}^{\ell}$ is the solution of

$\left(\left(\mathcal{A}_{x}^{l}\mathcal{A}_{z}^{k}\left(K(x,z)\right)\right)_{l,k=1,\dots,\ell}+\lambda\mathcal{I}_{\mathbb{R}^{\ell\times\ell}}\right)\alpha=g,$
(2.15)

then the scalar spherical spline function corresponding to the coefficient
vector $\alpha$ is the unique minimizer of the corresponding Tikhonov
functional, that is

$S=\operatorname*{arg\,min}_{F\in{\mathscr{H}}}\left(\mathchoice{\left\lVert
g-\mathcal{A}F\right\rVert_{2}^{2}}{\lVert
g-\mathcal{A}F\rVert_{2}^{2}}{\lVert g-\mathcal{A}F\rVert_{2}^{2}}{\lVert
g-\mathcal{A}F\rVert_{2}^{2}}+\lambda\mathchoice{\left\lVert
F\right\rVert_{\mathscr{H}}^{2}}{\lVert F\rVert_{\mathscr{H}}^{2}}{\lVert
F\rVert_{\mathscr{H}}^{2}}{\lVert F\rVert_{\mathscr{H}}^{2}}\right).$

Note that $\mathcal{I}_{\mathbb{R}^{\ell\times\ell}}$ stands for the
$\ell\times\ell$-identity matrix.

##  3\. Functional Vector-valued Reproducing Kernel Hilbert Space Splines on
the Ball

In the previous section, we recapitulated that the scalar reproducing kernel
based spline approximation has several advantages. On the one hand, we have
the interpolation and best-approximation properties. On the other hand, by
means of the spline approximation the problem of minimizing the regularized
Tikhonov-Philips functional over an infinite-dimensional Hilbert space reduces
to solving a finite dimensional system of linear equations without loss of
information. We want to conserve these two outstanding properties for our
vector-valued RKHS splines. In order to do so, we start with the construction
of vector-valued Sobolev spaces over the ball, which is based on vector
Sobolev spaces on the sphere, [25], and scalar ones on the ball, [38, 1, 40].

###  3.1. Vector Sobolev Spaces on the Ball

As in the scalar-valued case, we want to construct vector splines based on
reproducing kernels over the ball. This method can easily be transferred to
arbitrary Hilbert spaces if a complete orthonormal system therein is known.

Based on [26, Ch. 12.4] and the idea of the spherical Helmholtz decomposition,
we split the space of all vector-valued arbitrarily often continuously
differentiable functions on the ball into three
$\mathchoice{{\left\langle\cdot,\cdot\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle\cdot,\cdot\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle\cdot,\cdot\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle\cdot,\cdot\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}$-orthogonal
subspaces. Due to knowledge of the null spaces of the operators related to the
functional inverse MEG and EEG problem, this decomposition reflects the
structure of the neuronal current. They are based on the spaces

$\mathrm{C}^{(i),\infty}(B_{R},\mathbb{R}^{3})\coloneqq\mathchoice{\left\\{f\in\mathrm{C}^{\infty}(B_{R},\mathbb{R}^{3})\;\middle|\;\tilde{O}^{(\iota)}_{\xi}f(x)=0\
\text{ if }\\{1,2,3\\}\ni\iota\neq
i\right\\}}{\left.\\{f\in\mathrm{C}^{\infty}(B_{R},\mathbb{R}^{3})\;\middle|\;\tilde{O}^{(\iota)}_{\xi}f(x)=0\
\text{ if }\\{1,2,3\\}\ni\iota\neq
i\\}\right.}{\left.\\{f\in\mathrm{C}^{\infty}(B_{R},\mathbb{R}^{3})\;\middle|\;\tilde{O}^{(\iota)}_{\xi}f(x)=0\
\text{ if }\\{1,2,3\\}\ni\iota\neq
i\\}\right.}{\left.\\{f\in\mathrm{C}^{\infty}(B_{R},\mathbb{R}^{3})\;\middle|\;\tilde{O}^{(\iota)}_{\xi}f(x)=0\
\text{ if }\\{1,2,3\\}\ni\iota\neq i\\}\right.},\qquad i\in\\{1,2,3\\},$

where $\tilde{O}^{(i)}$ is the adjoint operator of $\tilde{o}^{(i)}$ (with
respect to the $\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})$-norm). However, for the
construction of vector-valued splines, it is only necessary that the three
directions are $\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})$-orthogonal. Now, let one
direction indicated by the superscript $i\in\\{1,2,3\\}$ be arbitrary but fix.
For each subspace, we construct a family of orthonormal basis function
$\\{h_{m,n,j}^{(i)}\\}_{m,n,j}$ which can be separated into a radial and an
angular part. Note that for our particular application the orthonormal basis
functions are given in Thm. 2.1. In addition, we assume that the estimate

$\sup_{x\in B_{R}}\sum_{j=1}^{2n+1}\mathchoice{\left\lvert
h_{m,n,j}^{(i)}(x)\right\rvert^{2}}{\lvert
h_{m,n,j}^{(i)}(x)\rvert^{2}}{\lvert h_{m,n,j}^{(i)}(x)\rvert^{2}}{\lvert
h_{m,n,j}^{(i)}(x)\rvert^{2}}\eqqcolon B_{m,n}^{(i)}<\infty$

holds true for every $i\in\\{1,2,3\\}$, $m\in\mathbb{N}_{0}$, and
$n\in\mathbb{N}_{0_{i}}$.

######  Definition 3.1 (Vector Sobolev Space).

Let $R>0$ be a given radius, $i\in\\{1,2,3\\}$ be fixed, and let
$a^{(i)}\coloneqq(a_{m,n}^{(i)})_{m,n}$ be a given real sequence. We define a
functional
$\mathscr{E}^{(i)}\colon\mathrm{C}^{(i),\infty}(B_{R},\mathbb{R}^{3})\to\mathbb{R}$
by

$\mathscr{E}^{(i)}(f)\coloneqq\sum_{(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}}\sum_{j=1}^{2n+1}\left(a_{m,n}^{(i)}\right)^{2}\mathchoice{{\left\langle
f,{h}_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}^{2}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}^{2}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}^{2}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}^{2}}}$

and, consequently, the space $\mathscr{E}^{(i)}(a^{(i)},B_{R})$ is given by

$\mathscr{E}^{(i)}(a^{(i)},B_{R})\coloneqq\mathchoice{\left\\{f\in\mathrm{C}^{(i),\infty}(B_{R},\mathbb{R}^{3})\;\middle|\;\mathscr{E}^{(i)}(f)<\infty\text{
and }\mathchoice{{\left\langle
f,{h}_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}=0\text{ if
}a_{m,n}^{(i)}=0\right\\}}{\left.\\{f\in\mathrm{C}^{(i),\infty}(B_{R},\mathbb{R}^{3})\;\middle|\;\mathscr{E}^{(i)}(f)<\infty\text{
and }\mathchoice{{\left\langle
f,{h}_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}=0\text{ if
}a_{m,n}^{(i)}=0\\}\right.}{\left.\\{f\in\mathrm{C}^{(i),\infty}(B_{R},\mathbb{R}^{3})\;\middle|\;\mathscr{E}^{(i)}(f)<\infty\text{
and }\mathchoice{{\left\langle
f,{h}_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}=0\text{ if
}a_{m,n}^{(i)}=0\\}\right.}{\left.\\{f\in\mathrm{C}^{(i),\infty}(B_{R},\mathbb{R}^{3})\;\middle|\;\mathscr{E}^{(i)}(f)<\infty\text{
and }\mathchoice{{\left\langle
f,{h}_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,{h}_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}=0\text{ if
}a_{m,n}^{(i)}=0\\}\right.}.$

The space is equipped with the inner product defined for all $f$,
$g\in\mathscr{E}^{(i)}(a^{(i)},B_{R})$ by

$\mathchoice{{\left\langle f,g\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle
f,g\rangle_{\mathscr{H}^{(i)}}}}{{\langle
f,g\rangle_{\mathscr{H}^{(i)}}}}{{\langle
f,g\rangle_{\mathscr{H}^{(i)}}}}\coloneqq\sum_{(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}}\sum_{j=1}^{2n+1}{\left(a_{m,n}^{(i)}\right)}^{2}\mathchoice{{\left\langle
f,h^{(i)}_{m,n,j}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h^{(i)}_{m,n,j}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h^{(i)}_{m,n,j}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h^{(i)}_{m,n,j}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}\mathchoice{{\left\langle
g,h^{(i)}_{m,n,j}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
g,h^{(i)}_{m,n,j}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
g,h^{(i)}_{m,n,j}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
g,h^{(i)}_{m,n,j}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}.$

We call the completion of $\mathscr{E}^{(i)}(a^{(i)},B_{R})$ with respect to
$\mathchoice{\left\lVert\cdot\right\rVert_{\mathscr{H}^{(i)}}}{\lVert\cdot\rVert_{\mathscr{H}^{(i)}}}{\lVert\cdot\rVert_{\mathscr{H}^{(i)}}}{\lVert\cdot\rVert_{\mathscr{H}^{(i)}}}$
the _Sobolev space_ $\mathscr{H}^{(i)}(a^{(i)},B_{R})$. In addition, we define
with $a\coloneqq(a^{(1)},a^{(2)},a^{(3)})$ the space

$\mathscr{H}\coloneqq\mathscr{H}(a,B_{R},\mathbb{R}^{3})\coloneqq\bigoplus_{i=1}^{3}\mathscr{H}^{(i)}(a^{(i)},B_{R}).$

As in the scalar case, the smoothness of the functions $f\in\mathscr{H}$ is
closely related to the behavior of the generating sequence $a$. The increasing
behavior of the sequence needs to be compensated by a proportionate decay of
the generalized Fourier coefficients of the functions. This has consequences
for the smoothness of the functions. Accordingly, there exists an inclusion of
two different Hilbert spaces.

######  Corollary 3.2.

Let $i\in\\{1,2,3\\}$ be arbitrary. Let
$a^{(i)}\coloneqq(a_{m,n}^{(i)})_{m,n}$ and
$b^{(i)}\coloneqq(b_{m,n}^{(i)})_{m,n}$ be two real sequences with
$|a_{m,n}^{(i)}|\leq|b_{m,n}^{(i)}|$ for all $m\in\mathbb{N}_{0}$ and
$n\in\mathbb{N}_{0_{i}}$. Then

$\mathscr{H}^{(i)}\left(b^{(i)},{B_{R}}\right)\subset\mathscr{H}^{(i)}\left(a^{(i)},B_{R}\right).$

######  Definition 3.3 (Summability).

Let $i\in\\{1,2,3\\}$ be arbitrary. Let
$a^{(i)}\coloneqq(a_{m,n}^{(i)})_{m,n}$ and
$b^{(i)}\coloneqq(b_{m,n}^{(i)})_{m,n}$ be two given real sequences. The
sequence $a^{(i)}$ is said to be _$b^{(i)}$ -summable_ if

$\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
a_{m,n}^{(i)}\neq
0\end{subarray}}\left(\frac{b_{m,n}^{(i)}}{a_{m,n}^{(i)}}\right)^{2}\left(B_{m,n}^{(i)}\right)^{2}<\infty.$

If each element of the sequence $b^{(i)}$ is equal to one, $a^{(i)}$ is said
to be _summable_.

Note that the summability condition is closely related to a bound for the
orthonormal basis functions. Thus, if another orthonormal system is used for
the construction of the vector Sobolev space, then this summability condition
needs to be adapted.

######  Theorem 3.4.

Let for each $i\in\\{1,2,3\\}$ the sequence $a^{(i)}$ be $b^{(i)}$-summable,
where $b_{m,n}^{(i)}\neq 0$ for all $m\in\mathbb{N}_{0}$ and
$n\in\mathbb{N}_{0_{i}}$. Then each function
$f\in\mathscr{H}(a/b,B_{R},\mathbb{R}^{3})$ is also continuous and has a
uniformly convergent expansions in the basis $\\{h_{m,n,j}^{(i)}\\}$.

In this context, the quotient $a^{(i)}/b^{(i)}$ is understood as the element-
wise division of the sequences $a^{(i)}$ and $b^{(i)}$, that is
$a^{(i)}/b^{(i)}\coloneqq(a_{m,n}^{(i)}/b_{m,n}^{(i)})_{m,n}$.

###### Proof of Thm. 3.4.

Each function $f\in\mathscr{H}(a/b,B_{R},\mathbb{R}^{3})$ can be represented
by a Fourier series converging with respect to the corresponding
$\mathchoice{\left\lVert\cdot\right\rVert_{\mathscr{H}}}{\lVert\cdot\rVert_{\mathscr{H}}}{\lVert\cdot\rVert_{\mathscr{H}}}{\lVert\cdot\rVert_{\mathscr{H}}}$-norm.
The next estimate proves that the iterated Fourier series also converges
uniformly, due to the Cauchy-Schwarz inequality for series and the bound of
the orthonormal basis functions. Thus, the following estimate holds true for
each $i\in\\{1,2,3\\}$ and all $x\in B_{R}$:

$\mathchoice{\left\lvert\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
m+n\geq M\end{subarray}}\sum_{j=1}^{2n+1}\mathchoice{{\left\langle
f,h_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}h_{m,n,j}^{(i)}(x)\right\rvert^{2}}{\lvert\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
m+n\geq M\end{subarray}}\sum_{j=1}^{2n+1}\mathchoice{{\left\langle
f,h_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}h_{m,n,j}^{(i)}(x)\rvert^{2}}{\lvert\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
m+n\geq M\end{subarray}}\sum_{j=1}^{2n+1}\mathchoice{{\left\langle
f,h_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}h_{m,n,j}^{(i)}(x)\rvert^{2}}{\lvert\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
m+n\geq M\end{subarray}}\sum_{j=1}^{2n+1}\mathchoice{{\left\langle
f,h_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}h_{m,n,j}^{(i)}(x)\rvert^{2}}\\\
\begin{aligned}
&\leq\left(\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
m+n\geq
M\end{subarray}}\sum_{j=1}^{2n+1}\left(\frac{a_{m,n}^{(i)}}{b_{m,n}^{(i)}}\right)^{2}\mathchoice{{\left\langle
f,h_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}^{2}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}^{2}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}^{2}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}^{2}}}\right)\left(\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
m+n\geq
M\end{subarray}}\left(\frac{b_{m,n}^{(i)}}{a_{m,n}^{(i)}}\right)^{2}\sum_{j=1}^{2n+1}\mathchoice{\left\lvert
h_{m,n,j}^{(i)}(x)\right\rvert^{2}}{\lvert
h_{m,n,j}^{(i)}(x)\rvert^{2}}{\lvert h_{m,n,j}^{(i)}(x)\rvert^{2}}{\lvert
h_{m,n,j}^{(i)}(x)\rvert^{2}}\right)\\\ &\leq\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}}^{2}}{\lVert f\rVert_{\mathscr{H}}^{2}}{\lVert
f\rVert_{\mathscr{H}}^{2}}{\lVert
f\rVert_{\mathscr{H}}^{2}}\left(\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
m+n\geq
M\end{subarray}}\left(\frac{b_{m,n}^{(i)}}{a_{m,n}^{(i)}}\right)^{2}\left(B_{m,n}^{(i)}\right)^{2}\right)\end{aligned}$

The right-hand side converges to zero as $M\to\infty$ due to the summability
condition. Hence, this iterated series converges uniformly and each summand is
continuous in $B_{R}$ as assumed. This estimate also implies the absolute
convergence of the iterated series. Due to Cauchy’s theorem of double series,
the absolute convergence of the iterated series suffices for the absolute
convergence of the double series and their limits coincide. ∎

###  3.2. Vector Splines on the Ball

For the construction of vector-valued reproducing kernel based splines, we
first need to construct appropriate reproducing kernels. For this method, they
are constructed via the tensor product of vector-valued orthonormal basis
functions. This approach has already been used in [26, 25] for vector-valued
splines on the unit sphere. Thus, the reproducing kernel
$\mathfrak{k}^{(i)}\colon B_{R}\times B_{R}\to\mathbb{R}^{3\times 3}$
considered in this section depends on $i\in\\{1,2,3\\}$ and is defined by

$\mathfrak{k}^{(i)}(x,y)\coloneqq\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
\kappa_{m,n}^{(i)}\neq
0\end{subarray}}\sum_{j=1}^{2n+1}\left(\kappa_{m,n}^{(i)}\right)^{-2}h_{m,n,j}^{(i)}(x)\otimes
h_{m,n,j}^{(i)}(y),\qquad x,\,y\in B_{R}.$ (3.1)

Reproducing kernels of this kind are also introduced in [6, 10] and used for
the construction of interpolating vector splines and vector-valued wavelets on
the sphere. The Sobolev space
$\mathscr{H}^{(i)}((\kappa_{m,n}^{(i)})_{m,n},B_{R},\mathbb{R}^{3})$ is the
natural Sobolev space containing the reproducing kernel with argument fixed,
that is
$\mathfrak{k}^{(i)}(x,\cdot)\in\mathscr{H}^{(i)}((\kappa_{m,n}^{(i)})_{m,n},B_{R},\mathbb{R}^{3})$
and
$\mathfrak{k}^{(i)}(\cdot,y)\in\mathscr{H}^{(i)}((\kappa_{m,n}^{(i)})_{m,n},B_{R},\mathbb{R}^{3})$
for all $x,\ y\in B_{R}$.

######  Lemma 3.5.

Let $i\in\\{1,2,3\\}$ be fixed. If the sequence $(\kappa_{m,n}^{(i)})_{m,n}$
satisfies the summability condition, Def. 3.3, the series in (3.1) converges
pointwise with respect to the $2$-norm of matrices.

###### Proof.

The $2$-norm of rank-one outer products of two vectors $u$,
$v\in\mathbb{R}^{3}$ fulfills the estimate $|u\otimes v|\coloneqq\|u\otimes
v\|_{2}\leq\|u\|_{2}\|v\|_{2}$. Translated in our notation we achieve for the
tensor product of the orthonormal basis functions

$\left|h_{m,n,j}^{(i)}(x)\otimes
h_{m,n,j}^{(i)}(y)\right|\leq\left|h_{m,n,j}^{(i)}(x)\right|\left|h_{m,n,j}^{(i)}(y)\right|\qquad\text{
for all }x,\ y\in B_{R}.$

Hence, for the reproducing kernel we have a convergent majorant which is given
by

$\displaystyle|\mathfrak{k}^{(i)}(x,y)|$
$\displaystyle\leq\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
\kappa_{m,n}^{(i)}\neq
0\end{subarray}}\left(\kappa_{m,n}^{(i)}\right)^{-2}\sum_{j=1}^{2n+1}|h_{m,n,j}^{(i)}(x)||h_{m,n,j}^{(i)}(y)|$
$\displaystyle\leq\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
\kappa_{m,n}^{(i)}\neq
0\end{subarray}}\left(\kappa_{m,n}^{(i)}\right)^{-2}\left(\sum_{j=1}^{2n+1}|h_{m,n,j}^{(i)}(x)|^{2}\right)^{1/2}\left(\sum_{j=1}^{2n+1}|h_{m,n,j}^{(i)}(y)|^{2}\right)^{1/2}$
$\displaystyle\leq\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
\kappa_{m,n}^{(i)}\neq
0\end{subarray}}\left(\kappa_{m,n}^{(i)}\right)^{-2}\sup_{x\in
B_{R}}\sum_{j=1}^{2n+1}|h_{m,n,j}^{(i)}(x)|^{2}\leq\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
\kappa_{m,n}^{(i)}\neq
0\end{subarray}}\left(\kappa_{m,n}^{(i)}\right)^{-2}\left(B_{m,n}^{(i)}\right)^{2}<\infty.$
∎

For more details on vector-valued reproducing kernels, such as the reproducing
property, we refer to [25] since the conversion from the spherical case to the
ball case is straightforward. Furthermore, we assume that the sequence
$(\kappa_{m,n}^{(i)})_{m,n}$ is given in such a way that
$\mathscr{H}^{(i)}\coloneqq\mathscr{H}^{(i)}((\kappa_{m,n}^{(i)})_{m,n},B_{R})\subset\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})$.

The spline constructed in [25, 6] is used for an interpolation problem. In
contrast, we want to use the spline for approximating the solution of a
(functional) inverse problem. Having the vector-valued reproducing kernels at
hand, we develop in the following a novel approach to vector-valued splines
for inverse problems.

For the construction of approximation splines for functional inverse problems,
the data $g=(g_{1},\dots,g_{\ell})\in\mathbb{R}^{\ell}$ needs to be of the
form

$g=\mathcal{A}f\qquad\Leftrightarrow\qquad g_{k}=\mathcal{A}^{k}f,\qquad
k=1,\dots,\ell$

with the linear and continuous operator
$\mathcal{A}\coloneqq(\mathcal{A}^{1},\dots,\mathcal{A}^{\ell})^{\mathrm{T}}$,
the functionals $\mathcal{A}^{k}\colon\mathscr{H}^{(i)}\to\mathbb{R}$ for all
$k=1,\dots,\ell$, and the (sought) quantity $f\in\mathscr{H}^{(i)}$. In [38,
Ch. 6.4., Ch. 10], this method is extended to a scalar spline approximation on
the ball. Thus, we call a function of the form

$s(x)=\sum_{k=1}^{\ell}\alpha_{k}\mathcal{A}_{z}^{k}\left(\mathfrak{k}^{(i)}(z,x)\right),\qquad
x\in B_{R}$ (3.2)

with the coefficients
$\alpha=(\alpha_{k})_{k=1,\dots,\ell}\in\mathbb{R}^{\ell}$ a _spline function
in $\mathscr{H}^{(i)}$ subject to $\mathcal{A}$_. The set of all these spline
functions is denoted by
$\mathrm{Spline}((\kappa_{m,n}^{(i)})_{m,n},\mathcal{A})$ which is an
$\ell$-dimensional space. For this definition to make sense, we need to define
what we understand by applying the functional $\mathcal{A}^{k}$ to a tensor
product of two vector-valued functions. For a tensorial function
$\mathfrak{t}(z,x)=\sum_{i,j=1}^{3}\mathfrak{t}_{i,j}(z,x)\varepsilon^{i}\otimes\varepsilon^{j}$
with scalar component functions $\mathfrak{t}_{i,j}$, we set

$\mathcal{A}_{z}^{k}\mathfrak{t}(z,x)\coloneqq\sum_{j=1}^{3}\mathcal{A}_{z}^{k}\left(\sum_{i=1}^{3}\mathfrak{t}_{i,j}(z,x)\varepsilon^{i}\right)\varepsilon^{j},$
(3.3)

where $\\{\varepsilon^{i}\\}_{i=1,2,3}$ is the standard basis of
$\mathbb{R}^{3}$. Eventually, the vector-valued spline function has the
representation

$s(x)=\sum_{k=1}^{\ell}\alpha_{k}\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
\kappa_{m,n}^{(i)}\neq
0\end{subarray}}\sum_{j=1}^{2n+1}\left(\kappa_{m,n}^{(i)}\right)^{-2}h_{m,n,j}^{(i)}(x)\mathcal{A}^{k}h_{m,n,j}^{(i)},\qquad
x\in B_{R}.$ (3.4)

The convergence of the stated series is implied by the summability condition,
the convergence of the reproducing kernel series, and the linearity and
continuity of the functionals.

Several useful properties of (scalar) splines (over the ball) have already
been known. For example, in [38, Thm. 10.13-14], two minimum properties of
scalar splines over the ball are proved. These statements also hold true in
the vector-valued case, which is proved for the particular setting of
$\mathscr{H}^{(i)}$ being the product space of two scalar Sobolev spaces in
[9, 7] using tensor-valued reproducing kernels. Now we adapt these statements
to our setting. For this purpose, three central properties need to be
verified.

######  Lemma 3.6.

Let the reproducing kernel be given as in (3.1), then

  1. (1)

a kind of reproducing property holds true for all $k=1,\dots,\ell$, that is

$\mathchoice{{\left\langle\mathcal{A}_{x}^{k}\left(\mathfrak{k}^{(i)}(x,\cdot)\right),f\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}_{x}^{k}\left(\mathfrak{k}^{(i)}(x,\cdot)\right),f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}_{x}^{k}\left(\mathfrak{k}^{(i)}(x,\cdot)\right),f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}_{x}^{k}\left(\mathfrak{k}^{(i)}(x,\cdot)\right),f\rangle_{\mathscr{H}^{(i)}}}}=\mathcal{A}^{k}f,$

for all $f\in\mathscr{H}^{(i)}$,

  2. (2)

for all $k,\,l=1,\dots,\ell$ the following relation holds true:

$\mathcal{A}^{l}_{x}\mathcal{A}^{k}_{z}\mathfrak{k}^{(i)}(x,z)=\mathchoice{{\left\langle\mathcal{A}^{k}_{x}\mathfrak{k}^{(i)}(\cdot,x),\mathcal{A}^{k}_{z}\mathfrak{k}^{(i)}(\cdot,z)\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{k}_{x}\mathfrak{k}^{(i)}(\cdot,x),\mathcal{A}^{k}_{z}\mathfrak{k}^{(i)}(\cdot,z)\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{k}_{x}\mathfrak{k}^{(i)}(\cdot,x),\mathcal{A}^{k}_{z}\mathfrak{k}^{(i)}(\cdot,z)\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{k}_{x}\mathfrak{k}^{(i)}(\cdot,x),\mathcal{A}^{k}_{z}\mathfrak{k}^{(i)}(\cdot,z)\rangle_{\mathscr{H}^{(i)}}}},$

and

  3. (3)

every spline function
$s\in\mathrm{Spline}((\kappa_{m,n}^{(i)})_{m,n},\mathcal{A})$ satisfies the
relation

$\mathchoice{{\left\langle s,f\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle
s,f\rangle_{\mathscr{H}^{(i)}}}}{{\langle
s,f\rangle_{\mathscr{H}^{(i)}}}}{{\langle
s,f\rangle_{\mathscr{H}^{(i)}}}}=\sum_{k=1}^{\ell}\alpha_{k}\mathcal{A}^{k}f$
(3.5)

for all $f\in\mathscr{H}^{(i)}$.

###### Proof.

We start with the proof of the first item. Then, we obtain with the definition
of the $\mathscr{H}^{(i)}$-inner product, (3.3), (3.4), and Thm. 3.4 as well
as the linearity and continuity of the $\mathcal{A}^{k}$

$\mathchoice{{\left\langle\mathcal{A}_{x}^{k}\left(\mathfrak{k}^{(i)}(x,\cdot)\right),f\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}_{x}^{k}\left(\mathfrak{k}^{(i)}(x,\cdot)\right),f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}_{x}^{k}\left(\mathfrak{k}^{(i)}(x,\cdot)\right),f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}_{x}^{k}\left(\mathfrak{k}^{(i)}(x,\cdot)\right),f\rangle_{\mathscr{H}^{(i)}}}}\begin{aligned}
&=\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
\kappa_{m,n}^{(i)}\neq
0\end{subarray}}\sum_{j=1}^{2n+1}\left(\kappa_{m,n}^{(i)}\right)^{2}\left(\kappa_{m,n}^{(i)}\right)^{-2}\mathcal{A}^{k}\left(h_{m,n,j}^{(i)}\right)\mathchoice{{\left\langle
f,h_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}\\\
&=\mathcal{A}^{k}\left(\sum_{\begin{subarray}{c}(m,n)\in\mathbb{N}_{0}\times\mathbb{N}_{0_{i}}\\\
\kappa_{m,n}^{(i)}\neq
0\end{subarray}}\sum_{j=1}^{2n+1}\mathchoice{{\left\langle
f,h_{m,n,j}^{(i)}\right\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}{{\langle
f,h_{m,n,j}^{(i)}\rangle_{\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})}}}h_{m,n,j}^{(i)}\right)\\\
&=\mathcal{A}^{k}f\end{aligned}$ (3.6)

The last step is valid since
$\mathscr{H}^{(i)}\subset\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})$ and, hence, $f$
can be represented by the Fourier series. The second item is a particular case
of the first with $f=\mathcal{A}^{k}_{z}\mathfrak{k}^{(i)}(\cdot,z)$.

For the proof of the last statement, we only need to use the representation of
the vector-valued spline from (3.2), the linearity of the inner-product, and
the first item of this lemma. ∎

With this preliminary work, the ideas of the original proofs, which can for
instance be found in [38], of the next statements are still valid in the
vector-valued case.

######  Theorem 3.7.

Let $g\in\mathbb{R}^{\ell}$ be the given data and the spline function
$s\in\mathscr{H}^{(i)}$ be unknown. Then the spline interpolation problem
$\mathcal{A}^{k}s=g_{k}$ for all $k=1,\dots,\ell$ is uniquely solvable if and
only if the functionals
$\mathchoice{\left\\{\mathcal{A}^{k}\right\\}_{k=1,\dots,\ell}}{\left.\\{\mathcal{A}^{k}\\}_{k=1,\dots,\ell}\right.}{\left.\\{\mathcal{A}^{k}\\}_{k=1,\dots,\ell}\right.}{\left.\\{\mathcal{A}^{k}\\}_{k=1,\dots,\ell}\right.}$
are linearly independent.

###### Proof.

Via (3.2), we see that the interpolation problem is equivalent to solving a
linear system of equations with the given data as the right-hand side and the
matrix

$\left(\mathcal{A}_{x}^{l}\mathcal{A}_{z}^{k}\left(\mathfrak{k}^{(i)}(z,x)\right)\right)_{l,k=1,\dots,\ell},$

which is uniquely solvable if and only if the matrix is regular. Via item (2)
of Lem. 3.6, we get that this matrix is a Gramian matrix. Item (1) of this
lemma provides us with the property that
$\mathchoice{\left\\{\mathcal{A}^{k}\right\\}_{k=1,\dots,\ell}}{\left.\\{\mathcal{A}^{k}\\}_{k=1,\dots,\ell}\right.}{\left.\\{\mathcal{A}^{k}\\}_{k=1,\dots,\ell}\right.}{\left.\\{\mathcal{A}^{k}\\}_{k=1,\dots,\ell}\right.}$
is linearly independent if and only if
$\mathchoice{\left\\{\mathcal{A}_{x}^{k}\mathfrak{k}^{(i)}(x,\cdot)\right\\}_{k=1,\dots,\ell}}{\left.\\{\mathcal{A}_{x}^{k}\mathfrak{k}^{(i)}(x,\cdot)\\}_{k=1,\dots,\ell}\right.}{\left.\\{\mathcal{A}_{x}^{k}\mathfrak{k}^{(i)}(x,\cdot)\\}_{k=1,\dots,\ell}\right.}{\left.\\{\mathcal{A}_{x}^{k}\mathfrak{k}^{(i)}(x,\cdot)\\}_{k=1,\dots,\ell}\right.}$
is linearly independent. ∎

Again with (3.5), the proofs of the following theorems are immediate
consequences of the ones in [38, Thm. 10.13-10.16] and therefore skipped for
the sake of brevity.

######  Theorem 3.8 (Minimum Properties).

Let $\mathscr{H}^{(i)}\subset\mathrm{L}^{2}(B_{R},\mathbb{R}^{3})$ be a given
Sobolev space and $\mathcal{A}^{k}\colon\mathscr{H}^{(i)}\to\mathbb{R}$ be
bounded linear functionals for all $k=1,\dots,\ell$ that are linearly
independent. Then the following properties hold true:

  1. (1)

If $g\in\mathbb{R}^{\ell}$ is a given vector and the spline $s$ is given by
$\mathcal{A}^{k}s=g_{k}$ for all $k=1,\dots,\ell$, then $s$ is the unique
minimizer of

$\mathchoice{\left\lVert s\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}}{\lVert s\rVert_{\mathscr{H}^{(i)}}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}}=\min\mathchoice{\left\\{\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}\;\middle|\;f\in\mathscr{H}^{(i)}\text{ with
}\mathcal{A}^{k}f=g_{k}\text{ for all
}k=1,\dots,\ell\right\\}}{\left.\\{\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}\;\middle|\;f\in\mathscr{H}^{(i)}\text{ with
}\mathcal{A}^{k}f=g_{k}\text{ for all
}k=1,\dots,\ell\\}\right.}{\left.\\{\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}\;\middle|\;f\in\mathscr{H}^{(i)}\text{ with
}\mathcal{A}^{k}f=g_{k}\text{ for all
}k=1,\dots,\ell\\}\right.}{\left.\\{\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}\;\middle|\;f\in\mathscr{H}^{(i)}\text{ with
}\mathcal{A}^{k}f=g_{k}\text{ for all }k=1,\dots,\ell\\}\right.}.$

  2. (2)

If $f\in\mathscr{H}^{(i)}$ is a given function and the spline $s$ is defined
by $\mathcal{A}^{k}s=\mathcal{A}^{k}f$ for all $k=1,\dots,\ell$, then $s$ is
the unique minimizer of

$\mathchoice{\left\lVert f-s\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-s\rVert_{\mathscr{H}^{(i)}}}{\lVert f-s\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-s\rVert_{\mathscr{H}^{(i)}}}=\min\mathchoice{\left\\{\mathchoice{\left\lVert
f-\bar{s}\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}\;\middle|\;\bar{s}\in\mathrm{Spline}\left((\kappa_{m,n}^{(i)})_{m,n},\mathcal{A}\right)\right\\}}{\left.\\{\mathchoice{\left\lVert
f-\bar{s}\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}\;\middle|\;\bar{s}\in\mathrm{Spline}\left((\kappa_{m,n}^{(i)})_{m,n},\mathcal{A}\right)\\}\right.}{\left.\\{\mathchoice{\left\lVert
f-\bar{s}\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}\;\middle|\;\bar{s}\in\mathrm{Spline}\left((\kappa_{m,n}^{(i)})_{m,n},\mathcal{A}\right)\\}\right.}{\left.\\{\mathchoice{\left\lVert
f-\bar{s}\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-\bar{s}\rVert_{\mathscr{H}^{(i)}}}\;\middle|\;\bar{s}\in\mathrm{Spline}\left((\kappa_{m,n}^{(i)})_{m,n},\mathcal{A}\right)\\}\right.}.$

######  Theorem 3.9 (Shannon Sampling Theorem).

Any spline function
$s\in\mathrm{Spline}\left((\kappa_{m,n}^{(i)})_{m,n},\mathcal{A}\right)$ is
representable by its samples $\mathcal{A}^{k}s$ as

$s(x)=\sum_{k=1}^{\ell}(\mathcal{A}^{k}s)L_{k}(x),\qquad x\in B_{R},$

where

$L_{k}(x)=\sum_{j=1}^{\ell}\alpha_{j}^{(k)}\mathcal{A}^{j}\mathfrak{k}^{(i)}(\cdot,x),\qquad
x\in B_{R},$

with $\alpha_{j}^{(k)}$ given as solutions of the linear equation systems

$\sum_{j=1}^{\ell}\alpha_{j}^{(k)}\mathcal{A}^{l}\mathcal{A}^{j}\mathfrak{k}^{(i)}(\cdot,\cdot)=\delta_{k,l}\qquad\text{
for all }k,\,l=1,\dots,N.$

######  Theorem 3.10 (Spline Approximation).

Let $g\in\mathbb{R}^{\ell}$ and a regularization parameter $\lambda>0$ be
given. Let the bounded linear functionals
$\mathcal{A}^{k}\colon\mathscr{H}^{(i)}\to\mathbb{R}$, with $k=1,\dots,\ell$,
be linearly independent. If the vector $\alpha=(\alpha_{k})_{k=1,\dots,\ell}$
is the solution of

$\left(\left(\mathcal{A}_{x}^{l}\mathcal{A}_{z}^{k}\left(\mathfrak{k}^{(i)}(z,x)\right)\right)_{l,k=1,\dots,\ell}+\lambda\mathcal{I}_{\mathbb{R}^{\ell\times\ell}}\right)\alpha=g,$
(3.7)

then the spline function given by (3.2) is the unique minimizer of the
corresponding Tikhonov-functional, that is

$s=\operatorname*{arg\,min}_{f\in\mathscr{H}^{(i)}}\left(\mathchoice{\left\lVert
g-\mathcal{A}f\right\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}+\lambda\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}^{2}}\right).$

###### Proof.

Let $s\in\mathrm{Spline}\left((\kappa_{m,n}^{(i)})_{m,n},\mathcal{A}\right)$
be the unique solution of (3.7) which exists since the occurring matrix is
positive definite, see Lem. 3.6. We first prove the auxiliary statement

$\mathchoice{{\left\langle\mathcal{A}^{\ast}g-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s,s-f\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}g-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s,s-f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}g-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s,s-f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}g-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s,s-f\rangle_{\mathscr{H}^{(i)}}}}=0$

for all $f\in\mathscr{H}^{(i)}$. For this purpose, we notice by means of the
spline definition, (3.2), and the canonical matrix-vector product that

$\mathcal{A}_{x}s(x)=\left(\mathcal{A}_{x}^{l}\mathcal{A}_{z}^{k}\left(\mathfrak{k}^{(i)}(z,x)\right)\right)_{l,k=1,\dots,\ell}\alpha.$

In addition, $\mathcal{A}^{\ast}$ can be applied to (3.7) and the resulting
equation can be equivalently reformulated to

$\mathcal{A}^{\ast}\mathcal{A}s+\lambda\mathcal{A}^{\ast}\alpha=\mathcal{A}^{\ast}g.$

Thus,

$\mathcal{A}^{\ast}g-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s=\mathcal{A}^{\ast}\mathcal{A}s+\lambda\mathcal{A}^{\ast}\alpha-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s=\lambda\left(\mathcal{A}^{\ast}\alpha-s\right).$

It remains to prove that the right-hand side of the last equation vanishes.
For all $f\in\mathscr{H}^{(i)}$ we have with item (3) of Lem. 3.6 the relation

$\mathchoice{{\left\langle\mathcal{A}^{\ast}\alpha-s,f\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}\alpha-s,f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}\alpha-s,f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}\alpha-s,f\rangle_{\mathscr{H}^{(i)}}}}=\mathchoice{{\left\langle\alpha,\mathcal{A}f\right\rangle_{\mathbb{R}^{\ell}}}}{{\langle\alpha,\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}{{\langle\alpha,\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}{{\langle\alpha,\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}-\mathchoice{{\left\langle
s,f\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle
s,f\rangle_{\mathscr{H}^{(i)}}}}{{\langle
s,f\rangle_{\mathscr{H}^{(i)}}}}{{\langle
s,f\rangle_{\mathscr{H}^{(i)}}}}=\mathchoice{{\left\langle\alpha,\mathcal{A}f\right\rangle_{\mathbb{R}^{\ell}}}}{{\langle\alpha,\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}{{\langle\alpha,\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}{{\langle\alpha,\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}-\alpha\cdot\mathcal{A}f=0.$

For the regularized Tikhonov-functional, we obtain for all
$f\in\mathscr{H}^{(i)}$ the estimate

$\displaystyle\phantom{=\ }\mathchoice{\left\lVert
g-\mathcal{A}f\right\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}+\lambda\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert f\rVert_{\mathscr{H}^{(i)}}^{2}}$
$\displaystyle=\mathchoice{\left\lVert
g-\mathcal{A}s+\mathcal{A}s-\mathcal{A}f\right\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s+\mathcal{A}s-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s+\mathcal{A}s-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s+\mathcal{A}s-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}+\lambda\mathchoice{\left\lVert
f-s+s\right\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f-s+s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f-s+s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f-s+s\rVert_{\mathscr{H}^{(i)}}^{2}}$ $\displaystyle=\mathchoice{\left\lVert
g-\mathcal{A}s\right\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}+\mathchoice{\left\lVert\mathcal{A}s-\mathcal{A}f\right\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert\mathcal{A}s-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert\mathcal{A}s-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert\mathcal{A}s-\mathcal{A}f\rVert_{\mathbb{R}^{\ell}}^{2}}+2\mathchoice{{\left\langle
g-\mathcal{A}s,\mathcal{A}s-\mathcal{A}f\right\rangle_{\mathbb{R}^{\ell}}}}{{\langle
g-\mathcal{A}s,\mathcal{A}s-\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}{{\langle
g-\mathcal{A}s,\mathcal{A}s-\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}{{\langle
g-\mathcal{A}s,\mathcal{A}s-\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}+\lambda\left(\mathchoice{\left\lVert
f-s\right\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f-s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f-s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
f-s\rVert_{\mathscr{H}^{(i)}}^{2}}+2\mathchoice{{\left\langle
f-s,s\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle
f-s,s\rangle_{\mathscr{H}^{(i)}}}}{{\langle
f-s,s\rangle_{\mathscr{H}^{(i)}}}}{{\langle
f-s,s\rangle_{\mathscr{H}^{(i)}}}}+\mathchoice{\left\lVert
s\right\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}\right)$
$\displaystyle\geq\mathchoice{\left\lVert
g-\mathcal{A}s\right\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}+\lambda\mathchoice{\left\lVert
s\right\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}+2\mathchoice{{\left\langle
g-\mathcal{A}s,\mathcal{A}s-\mathcal{A}f\right\rangle_{\mathbb{R}^{\ell}}}}{{\langle
g-\mathcal{A}s,\mathcal{A}s-\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}{{\langle
g-\mathcal{A}s,\mathcal{A}s-\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}{{\langle
g-\mathcal{A}s,\mathcal{A}s-\mathcal{A}f\rangle_{\mathbb{R}^{\ell}}}}+2\lambda\mathchoice{{\left\langle
f-s,s\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle
f-s,s\rangle_{\mathscr{H}^{(i)}}}}{{\langle
f-s,s\rangle_{\mathscr{H}^{(i)}}}}{{\langle
f-s,s\rangle_{\mathscr{H}^{(i)}}}}$ $\displaystyle=\mathchoice{\left\lVert
g-\mathcal{A}s\right\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}+\lambda\mathchoice{\left\lVert
s\right\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}+2\left(\mathchoice{{\left\langle\mathcal{A}^{\ast}(g-\mathcal{A}s),s-f\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}(g-\mathcal{A}s),s-f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}(g-\mathcal{A}s),s-f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}(g-\mathcal{A}s),s-f\rangle_{\mathscr{H}^{(i)}}}}-\mathchoice{{\left\langle\lambda
s,s-f\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle\lambda
s,s-f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\lambda
s,s-f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\lambda
s,s-f\rangle_{\mathscr{H}^{(i)}}}}\right)$
$\displaystyle=\mathchoice{\left\lVert
g-\mathcal{A}s\right\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}+\lambda\mathchoice{\left\lVert
s\right\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}+2\mathchoice{{\left\langle\mathcal{A}^{\ast}g-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s,s-f\right\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}g-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s,s-f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}g-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s,s-f\rangle_{\mathscr{H}^{(i)}}}}{{\langle\mathcal{A}^{\ast}g-\mathcal{A}^{\ast}\mathcal{A}s-\lambda
s,s-f\rangle_{\mathscr{H}^{(i)}}}}$ $\displaystyle=\mathchoice{\left\lVert
g-\mathcal{A}s\right\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}{\lVert
g-\mathcal{A}s\rVert_{\mathbb{R}^{\ell}}^{2}}+\lambda\mathchoice{\left\lVert
s\right\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert
s\rVert_{\mathscr{H}^{(i)}}^{2}}{\lVert s\rVert_{\mathscr{H}^{(i)}}^{2}}.$ ∎

Note that equality only holds if $s=f$.

######  Theorem 3.11 (Convergence Result).

Let $\mathchoice{\left\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\right\\}}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}$ be a countable, infinite, and linearly independent
system of linear and continuous functionals. For a fixed
$f\in\mathscr{H}^{(i)}$ and every $N\in\mathbb{N}$, let the spline function
$s_{N}$ be given by

$s_{N}\in\mathrm{Spline}\left(\left(\kappa_{m,n}^{(i)}\right)_{m,n},\left(\mathcal{A}^{1},\dots,\mathcal{A}^{N}\right)\right)\qquad\text{
with }\mathcal{A}^{k}s_{N}=\mathcal{A}^{k}f\quad\text{ for all }k=1,\dots,N.$

If $\operatorname{span}\mathchoice{\left\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\right\\}}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}$ is dense in the dual space
$(\mathscr{H}^{(i)})^{\ast}$ of $\mathscr{H}^{(i)}$ then

$\lim_{N\to\infty}\mathchoice{\left\lVert
f-s_{N}\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert
f-s_{N}\rVert_{\mathscr{H}^{(i)}}}=0.$

###### Proof.

We prove the convergence result in two steps. First, we prove the weak
convergence of the splines. Second, the convergence of
$\mathchoice{\left\lVert s_{N}\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert s_{N}\rVert_{\mathscr{H}^{(i)}}}$ to
$\mathchoice{\left\lVert f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}$ is proven. Combined we obtain the desired
result.

For the weak convergence we choose an arbitrary element
$\mathcal{T}\in(\mathscr{H}^{(i)})^{\ast}$. Since
$\operatorname{span}\mathchoice{\left\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\right\\}}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}$ is dense in the dual space, we find for all
$\varepsilon>0$ an $\tilde{N}\in\mathbb{N}$ and coefficients
$(\alpha_{k})_{k}$ such that

$\tilde{\mathcal{T}}\coloneqq\sum_{k=1}^{\tilde{N}}\alpha_{k}\mathcal{A}^{k}$

satisfies
$\mathchoice{\left\lVert\mathcal{T}-\tilde{\mathcal{T}}\right\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\mathcal{T}-\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\mathcal{T}-\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\mathcal{T}-\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}\leq\varepsilon$.
In addition, we obtain with the interpolation condition of the spline function
for all $N\geq\tilde{N}$ the relation

$\tilde{\mathcal{T}}f=\sum_{k=1}^{\tilde{N}}\alpha_{k}\mathcal{A}^{k}f=\sum_{k=1}^{\tilde{N}}\alpha_{k}\mathcal{A}^{k}s_{N}=\tilde{\mathcal{T}}s_{N}.$

This yields together with the first minimum norm property of Thm. 3.8

$\displaystyle\mathchoice{\left\lvert\mathcal{T}f-\mathcal{T}s_{N}\right\rvert}{\lvert\mathcal{T}f-\mathcal{T}s_{N}\rvert}{\lvert\mathcal{T}f-\mathcal{T}s_{N}\rvert}{\lvert\mathcal{T}f-\mathcal{T}s_{N}\rvert}$
$\displaystyle=\mathchoice{\left\lvert\mathcal{T}f-\tilde{\mathcal{T}}f+\tilde{\mathcal{T}}{s_{N}}-\mathcal{T}s_{N}\right\rvert}{\lvert\mathcal{T}f-\tilde{\mathcal{T}}f+\tilde{\mathcal{T}}{s_{N}}-\mathcal{T}s_{N}\rvert}{\lvert\mathcal{T}f-\tilde{\mathcal{T}}f+\tilde{\mathcal{T}}{s_{N}}-\mathcal{T}s_{N}\rvert}{\lvert\mathcal{T}f-\tilde{\mathcal{T}}f+\tilde{\mathcal{T}}{s_{N}}-\mathcal{T}s_{N}\rvert}$
$\displaystyle\leq\mathchoice{\left\lVert\mathcal{T}-\tilde{\mathcal{T}}\right\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\mathcal{T}-\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\mathcal{T}-\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\mathcal{T}-\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}(\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}+\mathchoice{\left\lVert
s_{N}\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert s_{N}\rVert_{\mathscr{H}^{(i)}}})$
$\displaystyle\leq 2\varepsilon\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}},$

which implies the weak convergence.

For the convergence of the norms, we use the representation, [52, p. 91],

$\mathchoice{\left\lVert f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}=\sup_{\begin{subarray}{c}\mathcal{T}\in\operatorname{span}\mathchoice{\left\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\right\\}}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}{\left.\\{\mathcal{A}^{k}\mid
k\in\mathbb{N}\\}\right.}\\\
\mathchoice{\left\lVert\mathcal{T}\right\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\mathcal{T}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\mathcal{T}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\mathcal{T}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}\leq
1\end{subarray}}\mathchoice{\left\lvert\mathcal{T}f\right\rvert}{\lvert\mathcal{T}f\rvert}{\lvert\mathcal{T}f\rvert}{\lvert\mathcal{T}f\rvert}.$

Thus, for all $\varepsilon>0$, we can find an $\tilde{N}\in\mathbb{N}$ and
coefficients $(\alpha_{k})_{k}$ such that the operator

$\tilde{\mathcal{T}}\coloneqq\sum_{k=1}^{\tilde{N}}\alpha_{k}\mathcal{A}^{k}\quad\text{
with
}\quad\mathchoice{\left\lVert\tilde{\mathcal{T}}\right\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}\leq
1$

satisfies

$\mathchoice{\left\lVert f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}\leq\mathchoice{\left\lvert\tilde{\mathcal{T}}f\right\rvert}{\lvert\tilde{\mathcal{T}}f\rvert}{\lvert\tilde{\mathcal{T}}f\rvert}{\lvert\tilde{\mathcal{T}}f\rvert}+\varepsilon.$

For all $N\geq\tilde{N}$ we obtain again with the first minimum norm property
the estimate

$\mathchoice{\left\lVert f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}-\varepsilon\leq\mathchoice{\left\lvert\tilde{\mathcal{T}}f\right\rvert}{\lvert\tilde{\mathcal{T}}f\rvert}{\lvert\tilde{\mathcal{T}}f\rvert}{\lvert\tilde{\mathcal{T}}f\rvert}=\mathchoice{\left\lvert\tilde{\mathcal{T}}s_{N}\right\rvert}{\lvert\tilde{\mathcal{T}}s_{N}\rvert}{\lvert\tilde{\mathcal{T}}s_{N}\rvert}{\lvert\tilde{\mathcal{T}}s_{N}\rvert}\leq\mathchoice{\left\lVert\tilde{\mathcal{T}}\right\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}{\lVert\tilde{\mathcal{T}}\rVert_{(\mathscr{H}^{(i)})^{\ast}}}\mathchoice{\left\lVert
s_{N}\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}\leq\mathchoice{\left\lVert
s_{N}\right\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}{\lVert
s_{N}\rVert_{\mathscr{H}^{(i)}}}\leq\mathchoice{\left\lVert
f\right\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}}{\lVert
f\rVert_{\mathscr{H}^{(i)}}}{\lVert f\rVert_{\mathscr{H}^{(i)}}},$

which implies the convergence of the norms. ∎

Now, we have the theoretical foundations at hand in order to solve the vector-
valued inverse MEG and EEG problems by means of RKHS splines.

##  4\. Foundations for Implementation

###  4.1. Scalar Splines for the Spherical MEG Problem

In [20], the presented scalar spline method is used for solving the scalar
inverse MEG problem. However, therein a layer density constraint is used to
achieve uniqueness of the inverse problem. Now, we want to use the minimum
norm condition instead. For this purpose, the method needs to be adapted.
Recall that the minimum norm condition is meant with respect to the entire
neuronal current and not with respect to a scalar-valued part of it. Since the
focus of this paper lies on the novel vector-valued spline method, we only
give a short summary of the formulae required for the scalar spline
implementation. We use the following set of functions for the radial part of
our orthonormal system

$G_{n}(r)\coloneqq\sqrt{\frac{2n+5}{\varrho_{0}^{3}}}\left(\frac{r}{\varrho_{0}}\right)^{n+1},\qquad
r\in[0,\varrho_{0}],\,n\in\mathbb{N}.$ (4.1)

Note that these functions fulfill the relation

$G_{n}^{\prime}(\varrho_{0})\varrho_{0}-(n-1)G_{n}(\varrho_{0})=2\sqrt{\frac{2n+5}{\varrho_{0}^{3}}},\qquad
n\in\mathbb{N}.$

In addition, due to the null space of our operator we choose
$\kappa_{0}\coloneqq 0$. Eventually, we can conclude that the linear
functionals stated in (2.13) applied to the corresponding reproducing kernel
yield

$\mathcal{A}^{k}_{\mathrm{m}}K(x,\cdot)=\frac{-\mu_{0}}{2\pi\varrho_{0}^{2}}\nu(y_{k})\cdot\left(\sum_{\begin{subarray}{c}n\in\mathbb{N}\\\
\kappa_{n}\neq
0\end{subarray}}\sqrt{\frac{(2n+1)(2n+5)^{2}}{n+1}}\kappa_{n}^{-2}\frac{r^{n+1}}{s_{k}^{n+2}}\tilde{p}_{n}^{(1)}(\eta_{k},\xi)\right).$
(4.2)

Each entry of the spline matrix is given for all $l$,
$k=1,\dots,{\ell_{\mathrm{M}}}$ by

$\mathcal{A}^{l}_{\mathrm{m},x}\mathcal{A}^{k}_{\mathrm{m},z}K(x,z)=\frac{\mu_{0}^{2}}{\varrho_{0}^{3}}\sum_{\begin{subarray}{c}n\in\mathbb{N}\\\
\kappa_{n}\neq
0\end{subarray}}\sum_{j=1}^{2n+1}\frac{4(2n+5)}{{(n+1)(2n+1)}\kappa_{n}^{2}}\left(\frac{\varrho_{0}^{2}}{s_{k}s_{l}}\right)^{n+2}\left(\nu(y_{k})\cdot\tilde{y}_{n,j}^{(1)}(\eta_{k})\right)\left(\nu(y_{l})\cdot\tilde{y}_{n,j}^{(1)}(\eta_{l})\right).$
(4.3)

An alternative representation, which is more suitable for implementation, is
stated in (B.1).

###  4.2. Vector Splines for the Spherical MEG and EEG problem

For the implementation of the vector spline method for the spherical MEG and
EEG problem, we have to determine the occurring functionals and reproducing
kernels with their corresponding vector Sobolev spaces.

We do not want to go into details with the precise vector-valued integral
representation of the forward operator and consider the SVD-based series
representation instead, derived in [35, 34]. Therein, it is verified that the
series representations converge absolutely and uniformly on the particular
domains. That is

$\displaystyle\mathcal{A}_{\mathrm{M}}^{k}J$
$\displaystyle=-\mu_{0}\sum_{n=1}^{\infty}\sum_{j=1}^{2n+1}\sqrt{\frac{n\varrho_{0}}{(2n+1)(2n+3)}}\mathchoice{{\left\langle
J,\tilde{g}^{(3)}_{0,n,j}(\varrho_{0};\cdot)\right\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}^{(3)}_{0,n,j}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}^{(3)}_{0,n,j}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}^{(3)}_{0,n,j}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}\left(\frac{\varrho_{0}}{s_{k}}\right)^{n+1}\frac{1}{s_{k}}\nu(y_{k})\cdot\tilde{y}_{n,j}^{(1)}(\eta_{k}),$
(4.4) $\displaystyle\mathcal{A}_{\mathrm{E}}^{k}J$
$\displaystyle=\sum_{n=1}^{\infty}\sum_{j=1}^{2n+1}\frac{1}{\sqrt{n\varrho_{0}}}\beta^{(L)}_{n}\mathchoice{{\left\langle
J,\tilde{g}_{0,n,j}^{(i)}(\varrho_{0};\cdot)\right\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}_{0,n,j}^{(i)}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}_{0,n,j}^{(i)}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}{{\langle
J,\tilde{g}_{0,n,j}^{(i)}(\varrho_{0};\cdot)\rangle_{\mathrm{L}^{2}(B_{\varrho_{0}},\mathbb{R}^{3})}}}\left({(n+1)}\left(\frac{s_{k}}{\varrho_{L}}\right)^{2n+1}+{n}\right)\left(\frac{\varrho_{0}}{s_{k}}\right)^{n+1}Y_{n,j}(\eta_{k}),$
(4.5)

where $y_{k}=s_{k}\eta_{k}$ for all $k=1,\dots,\ell_{\mathrm{M}}$ or
$k=1,\dots,\ell_{\mathrm{E}}$, respectively, denotes the sensor positions of
the measurement devices. Recall that a definition of the used orthonormal
basis function is stated in Thm. 2.1. Due to properties of the orthonormal
basis functions, we are able to calculate for all $m\in\mathbb{N}_{0},\
n\in\mathbb{N}$ the value of $B_{m,n}^{(i)}$, see (A.1) for more steps,

$B_{m,n}^{(i)}=\sup_{x\in
B_{\varrho_{L}}}\sum_{j=1}^{2n+1}\mathchoice{\left\lvert
h_{m,n,j}^{(i)}\right\rvert}{\lvert h_{m,n,j}^{(i)}\rvert}{\lvert
h_{m,n,j}^{(i)}\rvert}{\lvert
h_{m,n,j}^{(i)}\rvert}=\frac{\left(4m+2t_{n}^{(i)}+3\right)(2n+1)}{4\pi\varrho_{L}^{3}}\binom{m+t_{n}^{(i)}+1/2}{m}^{2}$
(4.6)

Now we need to determine appropriate RKHS with reproducing kernels. For this
purpose, we set the sequences
$\kappa^{(2)}\coloneqq(\kappa_{m,n}^{(2)})_{m,n}$ and
$\kappa^{(3)}\coloneqq(\kappa_{m,n}^{(3)})_{m,n}$ by

$\displaystyle\kappa_{m,n}^{(2)}$
$\displaystyle\coloneqq\kappa_{n}^{(2)}\delta_{m,0},$
$\displaystyle\kappa_{m,n}^{(3)}$
$\displaystyle\coloneqq\kappa_{n}^{(3)}\delta_{m,0}\qquad\text{for all
}n\in\mathbb{N},$

where we furthermore assume that $\kappa_{n}^{(i)}\neq 0$ and
$\kappa_{0}^{(i)}=0$ for all $i\in\\{2,3\\}$ and all $n\in\mathbb{N}$. This
reflects the structure of the integral operator null spaces. We define
accordingly the tensor-valued kernels

$\mathfrak{k}^{(i)}(x,y)\coloneqq\sum_{n=1}^{\infty}\sum_{j=1}^{2n+1}\left(\kappa_{n}^{(i)}\right)^{-2}\tilde{g}_{0,n,j}^{(i)}(\varrho_{0},x)\otimes\tilde{g}_{0,n,j}^{(i)}(\varrho_{0},y),\qquad
i\in\\{2,3\\}.$

Due to the choice of the sequence and the used orthonormal basis function, it
is guaranteed that the corresponding spline function does not violate the
minimum-norm condition of the neuronal current.

The summability conditions are given by

$\displaystyle\sum_{n=1}^{\infty}\left(\kappa_{n}^{(3)}\right)^{-2}\frac{n}{2n+3}\left(\frac{\varrho_{0}}{\varrho_{L}}\right)^{2n+2}$
$\displaystyle<\infty,$
$\displaystyle\sum_{n=1}^{\infty}\left(\kappa_{n}^{(2)}\right)^{-2}\frac{(2n+1)^{2}}{n}\left(\frac{\varrho_{0}}{\varrho_{L-1}}\right)^{2n+2}$
$\displaystyle<\infty.$

According to Def. 3.3, they consist of the RKHS sequence $\kappa_{n}^{(i)}$,
the bound of the orthonormal basis functions $B_{0,n}^{(i)}$, and a sequence
$b_{n}^{(i)}$ which is closely related to the singular values of the
operators, [34, Thm. 6.1, Thm. 6.3]. In order to use the presented vector-
valued spline method, we need to verify that the functionals of (4.4) and
(4.5) are continuous linear functionals mapping from
$\mathscr{H}^{(i)}(\kappa^{(i)},B_{\varrho_{0}},\mathbb{R}^{3})$ to
$\mathbb{R}$ for $i\in\\{2,3\\}$.

For proving the boundedness of the functionals, we mainly use the summability
conditions.

######  Theorem 4.1 (Boundedness of the Vector Spline Functionals).

Let the functionals mapping the vector-valued current $J$ onto the measured
quantities of the MEG and EEG device, be given by (4.4) and (4.5). Let the
Hilbert spaces
$\mathscr{H}^{(i)}(\kappa^{(i)},B_{\varrho_{0}},\mathbb{R}^{3})$ to
$\mathbb{R}$ for $i\in\\{2,3\\}$ of vector valued functions on the ball be
given. Then
$\mathcal{A}_{\mathrm{M}}^{k}\colon\mathscr{H}^{(3)}(\kappa^{(3)},B_{\varrho_{0}},\mathbb{R}^{3})\to\mathbb{R}$
and
$\mathcal{A}_{\mathrm{E}}^{k}\colon\mathscr{H}^{(2)}(\kappa^{(2)},B_{\varrho_{0}},\mathbb{R}^{3})\to\mathbb{R}$
are linear and continuous functionals.

###### Proof.

The proof can be found in Appendix B. ∎

Now we can determine the following expression required for the representation
of the corresponding MEG and EEG spline.

######  Lemma 4.2.

The MEG and EEG functionals applied to the tensor-valued reproducing kernels,
have the representations

$\displaystyle\mathcal{A}_{\mathrm{M}}^{k}\mathfrak{k}^{(3)}(\cdot,x)$
$\displaystyle=-\frac{\mu_{0}}{4\pi}\sum_{n=1}^{\infty}\left(\kappa_{n}^{(3)}\right)^{-2}\frac{r^{n}}{s_{k}^{n+2}}\left((\nu(y_{k})\cdot\eta_{k})(\xi\wedge\eta_{k})P_{n}^{\prime}(\xi\cdot\eta_{k})\phantom{\frac{1}{{n+1}}}\right.$
$\displaystyle\times\left.-\frac{1}{{n+1}}\left((\xi\wedge\eta_{k})P_{n}^{\prime\prime}(\xi\cdot\eta_{k})\left(\xi\cdot(\nu(y_{k})-(\nu(y_{k})\cdot\eta_{k})\eta_{k})\right)+P_{n}^{\prime}(\xi\cdot\eta_{k})\xi\wedge(\nu(y_{k})-(\nu(y_{k})\cdot\eta_{k})\eta_{k})\right)\right)$
$\displaystyle\mathcal{A}_{\mathrm{E}}^{k}\mathfrak{k}^{(2)}(\cdot,x)$
$\displaystyle=\frac{1}{4\pi}\sum_{n=1}^{\infty}\left(\kappa_{n}^{(2)}\right)^{-2}\frac{(2n+1)^{3/2}}{\sqrt{n}}\frac{r^{n-1}}{s_{k}^{n+1}}\beta^{(L)}_{n}\left({(n+1)}\left(\frac{s_{k}}{\varrho_{L}}\right)^{2n+1}+{n}\right)\tilde{p}_{n}^{(2)}(\xi;\eta_{k}).$

###### Proof.

The proof can be found in Appendix B. ∎

######  Theorem 4.3.

The entries of the vector spline matrices for the MEG and EEG case are given
by

$\displaystyle\mathcal{A}_{\mathrm{M}}^{l}\mathcal{A}_{\mathrm{M}}^{k}\left(\mathfrak{k}^{(3)}(\cdot,\cdot)\right)$
$\displaystyle=\frac{\mu_{0}^{2}}{\varrho_{0}}\sum_{n=1}^{\infty}\sum_{j=1}^{2n+1}{\frac{n\,(\kappa_{n}^{(3)})^{-2}}{(2n+1)(2n+3)}}\left(\frac{\varrho_{0}^{2}}{s_{l}s_{k}}\right)^{n+2}\left(\nu(y_{l})\cdot\tilde{y}_{n,j}^{(1)}(\eta_{l})\right)\left(\nu(y_{k})\cdot\tilde{y}_{n,j}^{(1)}(\eta_{k})\right),$
$\displaystyle\mathcal{A}_{\mathrm{E}}^{l}\mathcal{A}_{\mathrm{E}}^{k}\left(\mathfrak{k}^{(2)}(\cdot,\cdot)\right)$
$\displaystyle=\frac{1}{4\pi}\sum_{n=1}^{\infty}\left(\kappa_{n}^{(2)}\right)^{-2}\frac{2n+1}{{n\varrho_{0}}}\left(\beta^{(L)}_{n}\right)^{2}\left(\frac{\varrho_{0}^{2}}{s_{l}s_{k}}\right)^{n+1}\left({(n+1)}\left(\frac{s_{l}}{\varrho_{L}}\right)^{2n+1}+{n}\right)$
$\displaystyle\phantom{=}\times\left({(n+1)}\left(\frac{s_{k}}{\varrho_{L}}\right)^{2n+1}+{n}\right)P_{n}(\eta_{l}\cdot\eta_{k}).$

###### Proof.

The proof can be found in Appendix B. ∎

##  5\. Numerics

In this section, we present the numerical results achieved by the scalar and
vector based RKHS spline method. First we discuss the MEG case in more detail.
As stated before, there exist several approaches to solve this inverse
problem. On the one hand, the scalar approach is discussed, where the neuronal
current $J$ is decomposed by means of the Helmholtz decomposition. This
approach results in a relation between the scalar quantity $A^{(1)}$ and the
measured data. Via the scalar spline method the quantity $A^{(1)}$ can be
reconstructed from the data. Afterwards, this result can be transferred to a
reconstruction of the vector-valued current $J$ by the expression, [34, Sec.
20.2.4.],

$J(x)=-\frac{4\pi}{\varrho_{0}^{2}}\sum_{\begin{subarray}{c}n\in\mathbb{N}\\\
\kappa_{n}\neq
0\end{subarray}}\sum_{j=1}^{2n+1}\frac{(2n+3)(2n+5)}{(n+1)\sqrt{n(2n+1)}}\kappa_{n}^{-2}\left(\sum_{k=1}^{\ell_{\mathrm{M}}}\alpha_{k}\frac{r^{n}}{s_{k}^{n+2}}\nu(y_{k})\cdot\tilde{y}_{n,j}^{(1)}(\eta_{k})\right)\tilde{y}_{n,j}^{(3)}(\xi).$
(5.1)

For the implementation, the addition theorem should be used in order to get
rid of the summation over the spherical harmonics of order $j$. The particular
calculation can be found in the appendix, (B.4). On the other hand, by means
of the SVD-based method a direct relation between the vector-valued quantity
$J$ and its measured effect is achieved. Using the vector spline method, $J$
can directly be reconstructed from the data. Before we compare these two
approaches, we need to introduce the setup used for the numerics.

For solving the inverse MEG and EEG problem, we use the three-shell model,
which is the multiple-shell model for $L=3$. In addition, the radii of the
shells and their conductivities are assumed to be

$\displaystyle\varrho_{0}$ $\displaystyle=$0.071\text{\,}\mathrm{m}$,$
$\displaystyle\varrho_{1}$ $\displaystyle=$0.072\text{\,}\mathrm{m}$,$
$\displaystyle\varrho_{2}$ $\displaystyle=$0.079\text{\,}\mathrm{m}$,$
$\displaystyle\varrho_{3}$ $\displaystyle=$0.085\text{\,}\mathrm{m}$,$ (5.2)
$\displaystyle\sigma_{0}$
$\displaystyle=$0.330\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$,$
$\displaystyle\sigma_{1}$
$\displaystyle=$1.000\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$,$
$\displaystyle\sigma_{2}$
$\displaystyle=$0.042\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$,$
$\displaystyle\sigma_{3}$
$\displaystyle=$0.330\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$.$

These values are going back to [20] and are based on considerations from [32].
In addition, we use a set of real setting of sensor positions collected during
a clinical trial of the CBU, Cambridge, UK. They are plotted in relation to a
ball modeling the cerebrum in Fig. 1. Therein, we can see the irregular
distribution of the measurement positions. In the MEG case, there is a huge
data gap in the facial area. In the EEG case, almost all sensor positions are
located in the upper hemisphere.

Figure 1. Position of the MEG (left) and EEG (right) sensors around the
cerebrum modeled by a ball with radius $\varrho_{0}$

For our numerical tests, we create a realistic but synthetic test current by
means of a linear combination of two splines. Then, the forward functional is
applied to this spline for generating the data.

In order to avoid the inverse crime, we consider the following aspects:

  * •

The splines and its corresponding RKHS used for generating the data must be
different than those used for the inversion. In the case of generating the
data, we use the sequence

$\kappa_{n,\text{data}}^{-2}\coloneqq\frac{h_{\text{data}}^{2n+2}}{(2n+5)(2n+1)},\qquad
n\in\mathbb{N}_{0},$ (5.3)

for building the reproducing kernels. For the free parameter, we choose
$h_{\text{data}}=0.8$. Due to the Legendre series, the corresponding series
representation of the reproducing kernel (see also Sec. 4.1) has a closed form
based on the Legendre series, [26, Eq. (3.2.32)] , that is

$K(x,z)=\frac{1}{4\pi\varrho_{0}^{3}}\sum_{n=0}^{\infty}h_{\text{data}}^{2n+2}\left(\frac{r|z|}{\varrho_{0}^{2}}\right)^{n+1}P_{n}\left(\xi\cdot\frac{z}{|z|}\right)=\frac{1}{4\pi\varrho_{0}^{4}}\frac{hr}{\left|\frac{h_{\text{data}}}{\varrho_{0}}x-\frac{\varrho_{0}}{h_{\text{data}}|z|^{3}}z\right|}.$

  * •

We do not calculate the forward functional by the SVD-based method stated in
(4.2). Instead, we consider the integral representation of the forward
function from (2.13). More precisely,

$\begin{multlined}\phantom{adfsdaf}\mathcal{A}^{l}_{\mathrm{m},x}\mathcal{A}^{k}_{\mathrm{m},z}K(x,z)=\\\
\phantom{adfsd}\mu_{0}^{2}\int_{B_{\varrho_{0}}\times
B_{\varrho_{0}}}\Delta_{x}\left(|x|\Delta_{z}\left(|z|K(x,z)\right)\right)\left(\nu(y_{l})\cdot\left(\nabla_{y}K_{\mathrm{m}}(z,y)\right|_{y=y_{l}}\right)\left(\nu(y_{k})\cdot\left(\nabla_{y}K_{\mathrm{m}}(x,y)\middle)\right|_{y=y_{k}}\right)\,\mathrm{d}(x,z).\end{multlined}\phantom{adfsdaf}\mathcal{A}^{l}_{\mathrm{m},x}\mathcal{A}^{k}_{\mathrm{m},z}K(x,z)=\\\
\phantom{adfsd}\mu_{0}^{2}\int_{B_{\varrho_{0}}\times
B_{\varrho_{0}}}\Delta_{x}\left(|x|\Delta_{z}\left(|z|K(x,z)\right)\right)\left(\nu(y_{l})\cdot\left(\nabla_{y}K_{\mathrm{m}}(z,y)\right|_{y=y_{l}}\right)\left(\nu(y_{k})\cdot\left(\nabla_{y}K_{\mathrm{m}}(x,y)\middle)\right|_{y=y_{k}}\right)\,\mathrm{d}(x,z).$

Note that we interchanged the Laplacian with the integration. This is valid,
due to the integrability and smoothness of the kernel $K$, see [5, Cor. 16.3.]
and [35, Thm. 7.13]. In addition, by means of the aid of Mathematica [51] we
are able to state the identity

$\phantom{adfsd}\Delta_{z}\left(|z|K(x,z)\right)=\frac{h}{2\pi\varrho_{0}^{2}}|x|\left(3\frac{\varrho_{0}^{2}}{h^{2}}-4x\cdot
z+\left(\frac{h|x||z|}{\varrho_{0}}\right)^{2}\right)\left(\frac{\varrho_{0}^{2}}{h^{2}}-2x\cdot
z+\left(\frac{h|x||z|}{\varrho_{0}}\right)^{2}\right)^{-3/2}.$ (5.4)

  * •

The Laplacian with respect to $x$ is calculated numerically by a finite
difference method of accuracy order $8$, [22, Tab. 1].

  * •

The closed representation of the integral kernels gradient is used, [35, Thm.
15.15],

$\displaystyle\phantom{adh}4\pi\nabla_{y}K_{\mathrm{m}}(x,y)\coloneqq$
$\displaystyle\nabla_{y}\sum_{k=1}^{\infty}\frac{r^{k}}{s^{k+1}(k+1)}P_{k}(\xi\cdot\eta)$
(5.5) $\displaystyle=$
$\displaystyle\begin{cases}\left(\frac{1}{sr((\xi\cdot\eta)^{2}-1)}\left[\left(\frac{s-x\cdot\eta}{|x-y|}-1\right)\xi-\left(\frac{\xi\cdot
y-r}{|x-y|}-\xi\cdot\eta\right)\eta\right]+\frac{1}{s^{2}}\eta\right),&|\xi\cdot\eta|\neq
1,\\\ \left(\frac{1}{s(s-r)}+\frac{1}{s^{2}}\right)\eta,&\xi\cdot\eta=1,\\\
\left(\frac{1}{s(s+r)}-\frac{1}{s^{2}}\right)\eta,&\xi\cdot\eta=-1.\\\
\end{cases}$ (5.6)

  * •

The integration is performed by a quasi-Monte Carlo method, [15], with
$4\text{\times}{10}^{6}$ integration points distributed via a quasi-random
Kronecker sequence $t_{n}=\\{n\gamma\\}$,
$n\in\\{1,2,\dots,$4\text{\times}{10}^{6}$\\}$, and
$\gamma=(1/\Phi_{d}^{i})_{i=1,\dots,d}$ based on the $d$-dimensional
generalization of the golden ratio $\Phi_{d}$ as the solution of
$x^{d+1}=x+1$, [47];

  * •

We put additive Gaussian noise with $1$, $2$, $5$, and $10$ percent strength
on the data.

On the one hand, these steps lead us to a numerical method for calculating the
entries of the spline matrix. We will refer to this method as the numerical
analysis method, since it contains numerical differentiation and integration.
On the other hand, we have the possibility to calculate the spline matrix
based on the vector spherical harmonic expansion which is related to an SVD of
the operator, that is (4.3) combined with the addition theorem stated in
(B.1). This ansatz will be referred to as the SVD-based approach and yields to
a series representation for each entry. These series of Legendre polynomials
can be summed up efficiently by means of Clenshaw’s algorithm. Based on
Clenshaw’s algorithm, we have implemented a highly vectorized and optimized
code using MATLAB R2021a, [36]. For determining the run time of $0.4285\pm
0.0069\text{\,}\mathrm{s}$ in the case of $500$ summands per series, we have
run this algorithm $1000$ times sequentially on a local machine. In contrast,
this approach requires comprehensive knowledge of the forward operator in
order to obtain the series representation. This theoretical effort is rewarded
with a stable and very fast numerical method. For the numerical analysis
method, theoretical effort is required in order to gain the closed
representations of the kernels. Afterwards the numerical differentiation can
be operated fast by using vectorized code. However, each spline matrix entry
requires a numerical integration over $B_{\varrho_{0}}\times B_{\varrho_{0}}$.
In order to achieve an acceptable accuracy, $4\text{\times}{10}^{6}$
integration points are required for each integration over the three-
dimensional ball $B_{\varrho_{0}}$. We also vectorized and parallelized the
code and operated it on a node of the OMNI cluster (University of Siegen, 64
CPU cores, AMD EPYC 7452 at
$2.35\text{\,}\mathrm{GHz}3.35\text{\,}\mathrm{GHz}$,
$256\text{\,}\mathrm{GB}$ memory) in about $18$ days. This slowness reveals
the major disadvantage of this method, which comes along with the curse of
dimensions. However, the numerical integration may also offer the possibility
to extend this method to more real-shaped brain geometries, which is not
pursued within the scope of this paper.

In Fig. 2, the spline matrices achieved by the SVD-based method is visualized.
The matrix is ill-conditioned with condition number
$6.4053\text{\times}{10}^{11}$ with respect to the $2$-norm. In addition, the
element-wise relative derivation of the spline matrix achieved by the two
approaches described above is plotted on a logarithmic scale. Even though
there are several matrix entries with a higher point-wise deviation, the
average of the relative deviation over all matrix elements is
$0.13\text{\,}\mathrm{\char 37\relax}$. The higher point-wise deviation is
related to entries closer to zero. Hence, we conclude that the SVD-based
method as well as the numerical analysis based method yield comparable results
with respect to the accuracy.

Recall that the data generated by the linear combination of two splines
consists of the linear combination of the corresponding two rows of the spline
matrix. For generating the data, we calculate the spline matrix by means of
the numerical analysis method described before. For the inversion of the data,
we use the SVD-based method due to its immense speed-up to generate the spline
matrix. The used setting is described in the next section.

For a better comparison between the results obtained from the scalar and the
vector spline method, we use the same data set for both inversions.

Figure 2. Matrix of the scalar MEG splines with $h_{\text{data}}=0.8$ and
sequence (5.3), calculated by means of the SVD, (4.3), (left) and its
logarithmic plot, entry-wise relative difference to the matrix obtained by the
numerical method (right)

###  5.1. Spline Based Inversion of the MEG Data

As stated before, a linear and continuous functional
$\mathcal{A}^{k}_{\mathrm{m}}$ connecting the measured data $g_{k}$,
$k=1,\dots,\ell_{\mathrm{M}}$, with the radial component $A^{(1)}$ from the
vector field $a$ of the neuronal current Helmholtz decomposition (2.12) is
known, (4.2). In order to solve the problem

$\mathcal{A}^{k}_{\mathrm{m}}A^{(1)}=g_{k},\qquad
k=1,\dots,\ell_{\mathrm{M}},$

by means of the reproduced kernel based spline methods, we need to choose the
symbols of the reproducing kernel. We set

$\kappa_{n}^{-2}\coloneqq\frac{h^{n}}{n},\qquad
n\in\\{1,2,\dots,200\\},\,h\in\\{0.85,0.9,0.95,0.99\\}.$ (5.7)

The RKHS spline can now be determined as in (2.14). A visualization of splines
with only one and two non-vanishing coefficients $\alpha_{k}$ can be found in
Fig. 3. Therein the localization of the spline in an area around its node can
be seen. When the free parameter $h$ tends closer to $1$, the reproducing
kernel becomes more localized and the amplitude rises.

A Scalar RKHS spline

B Exact solution of spline based scalar test case for $A^{(1)}$

C Exact solution of ONB based scalar test case for $A^{(1)}$

Figure 3. Radially symmetric RKHS spline generated by the sequence (5.7) with
$h=0.95$ (left), a linear combination of two splines with $h=0.8$ and sequence
(5.3), which serves as the solution of $A^{(1)}$ in the spline based synthetic
test case (middle), and the ONB based test case solution $A^{(1)}$ (right).
All plotted on a sphere inside the cerebrum with radius $0.99\varrho_{0}$.

For each data set, we solve the regularized linear equation system, (2.15),
for $500$ different regularization parameters $\lambda$. The regularization
parameters are logarithmically uniformly distributed and weighted with the
absolute maximal value of the spline matrix. Hence, they are in the rage
$1\text{\times}{10}^{-25}1\text{\times}{10}^{-9}$. In order to pick the best
parameter, five parameter choice methods are applied afterwards, which are the
automatic and the manual L-curve method (LCM), the discrepancy principle (DP),
the quasi-optimality criterion (QOC), and the generalized cross validation
(GCV), see [29, 4, 3] and the references therein. For the resulting up to five
‘best’ parameters, the normalized root mean square error (NRMSE) of the
difference to the exact solution is calculated on the plotting grid. Then the
parameter choice method resulting in the smallest NRMSE is chosen. The L-curve
for the spline based synthetic test case with $5\text{\,}\mathrm{\char
37\relax}$ noise is plotted in Fig. 4.

Figure 4. Double logarithmic plot of the relative residual against the norm of
approximation depending on the regularization parameter. Regularization
parameters preferred by the five parameter choice methods are highlighted.

For each of the five different noise levels ($0\text{\,}\mathrm{\char
37\relax}$, $1\text{\,}\mathrm{\char 37\relax}$, $2\text{\,}\mathrm{\char
37\relax}$, $5\text{\,}\mathrm{\char 37\relax}$, $10\text{\,}\mathrm{\char
37\relax}$), we calculated the spline matrix with four different parameters
$h$, see (5.7). Among these $20$ inversions, the best regularization parameter
came only once from the GCV. In the other cases, the L-curve method provided
the best regularization parameter. With low noise levels (up to
$2\text{\,}\mathrm{\char 37\relax}$) the automatic L-curve method yields good
results. With higher noise levels, the a-priori range of the regularization
parameters gets more important. If the range is chosen to be wide, the system
of linear equations becomes easily over-regularized and the norm of the
solution is so small (magnitude of $1\text{\times}{10}^{-12}$ or less) such
that the automatic L-curve chooses regularization parameters leading to a
relative data misfit of $1\text{\times}{10}^{2}$ or higher. This behavior can
be avoided by restricting the parameter range a-priori or using the manual
L-curve method.

noise level | 0% | 1% | 2% | 5% | 10%
---|---|---|---|---|---
$h$ | $0.8500$ | $0.8500$ | $0.9900$ | $0.8500$ | $0.8500$
rel. residual | $4.2642\text{\times}{10}^{-08}$ | $0.0081$ | $0.0169$ | $0.0455$ | $0.0834$
rel. NRMSE $A^{(1)}$ | $5.3310\text{\times}{10}^{-04}$ | $0.0194$ | $0.0328$ | $0.0472$ | $0.0599$
rel. NRMSE $J$ | $5.3411\text{\times}{10}^{-04}$ | $0.0192$ | $0.0371$ | $0.0474$ | $0.0591$
$\lambda$ | $4.4054\text{\times}{10}^{-17}$ | $1.3395\text{\times}{10}^{-11}$ | $4.2240\text{\times}{10}^{-11}$ | $5.2497\text{\times}{10}^{-11}$ | $5.7572\text{\times}{10}^{-11}$
pcm | LCM man | LCM auto | LCM auto | LCM man | LCM man
$h$ | $0.3771$ | $0.3771$ | $0.3771$ | $0.3771$ | $0.3771$
rel. residual | $4.5957\text{\times}{10}^{-04}$ | $0.0080$ | $0.0166$ | $0.0451$ | $0.0874$
rel. NRMSE $J$ | $0.0051$ | $0.0140$ | $0.0240$ | $0.0416$ | $0.0455$
$\lambda$ | $3.6990\text{\times}{10}^{-15}$ | $1.4782\text{\times}{10}^{-15}$ | $1.4782\text{\times}{10}^{-15}$ | $2.5292\text{\times}{10}^{-14}$ | $1.4782\text{\times}{10}^{-14}$
pcm | LCM auto | LCM auto | LCM auto | LCM auto | LCM man
Table 1. Numerical results for the MEG synthetic spline based test case
achieved by the scalar spline method (top) and the vector spline method
(bottom)

Numerical results for the synthetic MEG data inversion such as the relative
residual as well as the relative NRMSE with respect to different noise levels
are listed in Tab. 1. For each noise level, we selected among the four values
of $h$, the one resulting in the best numerical result with respect to the
relative NRMSE. We can also see, that $h=0.85$ is chosen most frequently. A
possible explanation is that the corresponding width of the spline location
fits the best to the distances of the measurement positions and is the closest
to the synthetic test case value $h_{\text{data}}$. Besides the chosen $h$,
the relative data misfit or residual, respectively, is stated in Tab. 1. Even
with increasing noise level, the relative residual is often more than slightly
below the noise level.

Taking the relative NRMSE and visualizations of the reconstructions of
$A^{(1)}$ into account, Fig. 5, one can see that the approximations are not
over-fitted. In contrast, the active regions are covered well by the
reconstruction. In the case of lower noise levels, the amplitudes of the
current are matching. Only in the case of higher noise levels they become a
little bit blurred and faded, which goes back to the required stronger
regularization. Since deviations to the exact solution are hard to see for
lower noise levels, we have a closer look at them in Fig. 6. Near the active
regions, the deviations are the largest, which was anticipated. In addition,
larger differences can be seen in the lower part of the ball. This behavior is
directly related to the irregular distribution of measurement points, see Fig.
1, since there is only few information of the lower part of the ball and in
the facial area. The artifactual pattern in the deviation plots leads back to
the structure of the splines. All in all, the reconstructions of $A^{(1)}$
yield good and stable results with respect to increasing noise level.

A $0\text{\,}\mathrm{\char 37\relax}$ noise on the data

B $5\text{\,}\mathrm{\char 37\relax}$ noise on the data

C $10\text{\,}\mathrm{\char 37\relax}$ noise on the data

Figure 5. Reconstruction of the neuronal current component $A^{(1)}$ in the
synthetic spline based scalar MEG test case (3B) from $0\text{\,}\mathrm{\char
37\relax}$ (left), $5\text{\,}\mathrm{\char 37\relax}$ (middle), and
$10\text{\,}\mathrm{\char 37\relax}$ (right) noisy data

A $0\text{\,}\mathrm{\char 37\relax}$ noise on the data

B $5\text{\,}\mathrm{\char 37\relax}$ noise on the data

C $10\text{\,}\mathrm{\char 37\relax}$ noise on the data

Figure 6. Absolute deviation of the neuronal current component $A^{(1)}$ in
the synthetic spline based scalar MEG test case from $0\text{\,}\mathrm{\char
37\relax}$ (left), $5\text{\,}\mathrm{\char 37\relax}$ (middle),
$10\text{\,}\mathrm{\char 37\relax}$ (right) noisy data

Before we discuss the corresponding vector-valued reconstruction, we verify
the accuracy and stability of the numerical method by means of the inversion
of synthetic data which is not generated by a spline method at all. For this
purpose, we use the exact solution $A^{(1)}(r\xi)=r\xi\mapsto
0.1G_{3}(r)Y_{3,6}(\xi)$, which we refer to as the ONB based test case
solution and is plotted in Fig. 3C. The angular part of the function is given
by the a spherical harmonic of degree $3$ and order $6$, whereas the formula
for the radial part is stated in (4.1). For this particular test case, we have
also chosen $h=0.85$, since it yields the best results throughout our
numerical tests. Since the non-noisy and low-noisy reconstructions visually
coincide with the exact solution, we have only plotted the deviations in Fig.
7. Recall that the absolute maximum values of the exact solutions are
approximately $9.5$.

The approximation recovers the active regions and the amplitude well, even
with $10\text{\,}\mathrm{\char 37\relax}$ noise on the data, the
reconstruction does not seem to be blurry or fade out. This is also reflected
by the relative residual which is significantly smaller than the noise level
despite resulting in an satisfying small relative NRMSE. The particular values
are stated in the captions of Fig. 7. Besides, there are artifacts based on
the spline structure which can be seen in the deviation plots. Unattached from
the noise level, the areas most difficult to reconstruct are located in the
data gap area. This is especially revealed in the non-noisy case taking into
account that the exact solution is rotationally symmetric. Besides these small
artifacts in the data gap area, the used spline method can handle the
irregularly distributed data situation very well.

A noise level: $0\text{\,}\mathrm{\char 37\relax}$
rel. residual: $1.0092\text{\times}{10}^{-06}$
rel. NRMSE: $0.0076$

B noise level: $5\text{\,}\mathrm{\char 37\relax}$
rel. residual: $0.0313$
rel. NRMSE: $0.0623$

C noise level: $10\text{\,}\mathrm{\char 37\relax}$
rel. residual: $0.0619$
rel. NRMSE: $0.0891$

Figure 7. Absolute deviation of the reconstruction of the neuronal current
component $A^{(1)}$ in the ONB based synthetic scalar MEG test case to the
exact solution $r\xi\mapsto 0.1G_{3}(r)Y_{3,6}(\xi)$ depending on the noise
levels. For a better comparability, the colorbars were chosen equally.

However, we are mainly interested in the approximation and visualization of
the neuronal current, hence in the vector-valued quantity. In order to achieve
this, we first transfer the scalar spline solution to the neuronal current by
means of (5.1). The exact solutions of the current $J$ in the spline based as
well as the ONB based test case are plotted in Fig. 8. Besides the absolute
values of the currents, which are visualized via surface plots, the direction
of the currents are depicted via arrows.

Figure 8. Exact solution of current $J$ for the spline based synthetic test
case (left) and the ONB based test case (right)

We start with a look at the ONB based synthetic test case. In the non-noisy
case, the quality of the vector-valued current reconstruction is totally
comparable to the quality of the scalar one, Fig. 9. Also in the low-noisy
cases we observe that the main deviations to the exact solution are located in
the data gap area and the remaining approximation fits the exact solution very
well with respect of localization and amplitude. In the case of
$10\text{\,}\mathrm{\char 37\relax}$ noise on the data additional
reconstruction errors occur besides the data gap and the reconstruction
appears to be muted. This is also manifested in the values of the relative
NRMSE. Compared to the situation of the scalar approximation, see Fig. 7, the
relative NRMSE increases due to the transformation.

A Deviation from non-noisy data
rel. NRMSE: $0.0154$

B Deviation from $5\text{\,}\mathrm{\char 37\relax}$ noisy data
rel. NRMSE: $0.0787$

C Reconstruction from $10\text{\,}\mathrm{\char 37\relax}$ noise

D Deviation from $10\text{\,}\mathrm{\char 37\relax}$ noisy data
rel. NRMSE: $0.1033$

Figure 9. Numerical results for the ONB based synthetic test case achieved
from the transformation of the scalar reconstruction from Fig. 7. For a better
comparability, the colorbars of the deviation plots were chosen equally.

The relative NRMSE of the transferred vector-valued approximation of $J$ for
the spline based synthetic test case is also stated in Tab. 1 (top). In this
particular numerical experiment the relative NRMSE of the scalar
reconstruction and the transferred one are quite similar. However, they are
outperformed by the relative NRMSE of the direct vector-valued reconstruction
of $J$ via the vector spline method, whereas the relative residual remains
significantly lower than the noise level. This gain of quality in the
reconstruction can also be seen in the deviation plots, Fig. 10. In the first
row, the inversion of data with noise level $2\text{\,}\mathrm{\char
37\relax}$ is performed and with $10\text{\,}\mathrm{\char 37\relax}$ noise on
the data in the second row. In the left column the differences of the
transferred scalar reconstruction to the exact solution are plotted. In
comparison to the reconstructions of the direct vector spline method plotted
in the middle column, we see significantly higher deviations spread over the
entire ball. In the case of the vector spline reconstruction from
$2\text{\,}\mathrm{\char 37\relax}$, the highest deviations by far are located
in the data gap area. In the right column, the vector spline reconstructions
are plotted, which are stable with increasing noise level. Due to the
construction of the vector reproducing kernels, the vector splines are
significantly more localized than the scalar ones for higher parameter $h$.
Based on our numerical investigation, the adaption $h=0.85^{6}$ for the free
parameter yielded reasonable results in order to avoid an overfitting of the
data.

A Difference of converted scalar spline reconstruction
($2\text{\,}\mathrm{\char 37\relax}$ noise)

B Difference of vector spline reconstruction ($2\text{\,}\mathrm{\char
37\relax}$ noise)

C Reconstruction via vector splines ($2\text{\,}\mathrm{\char 37\relax}$
noise)

D Difference of converted scalar spline reconstruction
($10\text{\,}\mathrm{\char 37\relax}$ noise)

E Difference of vector spline reconstruction ($10\text{\,}\mathrm{\char
37\relax}$ noise)

F Reconstruction via vector splines ($10\text{\,}\mathrm{\char 37\relax}$
noise)

Figure 10. Deviation of the exact current to the converted reconstruction of
$J$ belonging to the scalar spline reconstruction of $A^{(1)}$ plotted in Fig.
3B (left column), to the direct reconstruction by means of the vector spline
method (middle column) and the reconstruction achieved by the vector splines
(right column) depending on $2\%$ noise (top row) and $10\%$ noise (bottom
row). For a better comparability, the colorbars in the left and middle column
were chosen equally.

We tested the vector spline method not only for the inverse MEG problem but
also for the inverse EEG problem. In this setting, we also generated synthetic
data by using a linear combination of two splines with $h_{\text{data}}=0.64$.
Having the inverse crime in mind, we used the sequence
$\kappa_{n,\text{data}}^{-2}\coloneqq nh^{n}$ for generating the data and
$\kappa_{n}^{-2}\coloneqq h^{n}$ with $h=0.85$ for the inversion. As an
additional obstacle, we located the test splines in the junction from the data
gap area to the area covered by sensors. The exact solution is plotted in Fig.
11A.

noise level | 0% | 1% | 2% | 5% | 10%
---|---|---|---|---|---
rel. residual | $3.0514\text{\times}{10}^{-09}$ | $0.0018$ | $0.0268$ | $0.0415$ | $0.0635$
rel. NRMSE $J$ | $0.0133$ | $0.0271$ | $0.0427$ | $0.0512$ | $0.0869$
$\lambda$ | $0.0233$ | $23.2964$ | $232.9644$ | $232.9644$ | $186.3716$
Table 2. Numerical results for the EEG synthetic test case achieved by the
vector spline method

Even though the synthetic setting is more difficult to solve as in the inverse
MEG test case, the results are comparably satisfying with respect to accuracy
and stability during increasing noise level, Tab. 2. Especially within small
noise levels, one can see in the deviation plots Fig. 11 (B) to (E) that the
reconstruction on the sensor covered hemisphere is more accurate than in the
data gap area, which was expected. For $10\text{\,}\mathrm{\char 37\relax}$
noise on the data this effect is smoothed away, since a higher regularization
was necessary to handle the noise. However, even for $10\text{\,}\mathrm{\char
37\relax}$ the reconstruction covers the active regions very satisfyingly and
is not over-fitting the data, see Fig. 11F .

A Exact spline based test current

B Deviation ($0\text{\,}\mathrm{\char 37\relax}$ noise)

C Deviation ($1\text{\,}\mathrm{\char 37\relax}$ noise)

D Deviation ($5\text{\,}\mathrm{\char 37\relax}$ noise)

E Deviation ($10\text{\,}\mathrm{\char 37\relax}$ noise)

F Reconstruction ($10\text{\,}\mathrm{\char 37\relax}$ noise)

Figure 11. Absolute deviations and reconstruction of neuronal current $J$ from
synthetic spline based EEG data with $0\text{\,}\mathrm{\char 37\relax}$,
$5\text{\,}\mathrm{\char 37\relax}$, and $10\text{\,}\mathrm{\char 37\relax}$
noise by using the vector spline method

Within our synthetic test cases for the inverse MEG and EEG problem, we have
seen that

  * •

the scalar spline method as well as the vector spline method is a stable
method for the solution of ill-posed inverse problems;

  * •

the spline methods are fast: the average of $700$ runs of the vector spline
method including parameter choice methods and plotting for the inverse MEG
problem yield $0.5654\pm 0.0261\text{\,}\mathrm{s}$ and $0.4285\pm
0.0069\text{\,}\mathrm{s}$ for the scalar spline method.

  * •

both methods are able to handle irregularly distributed and noisy data;

  * •

both methods keep the relative residual significantly below the noise level
without over-fitting the reconstruction;

  * •

the reconstructions coincides well with the exact solution even in the
presence of higher noise levels;

  * •

most deviations are located in the data gap area which was expected due to the
lack of information in this area;

  * •

the active regions are covered well and the amplitudes are reproduced
accurately in the non-noisy and low-noisy case (up to $5\text{\,}\mathrm{\char
37\relax}$) and very satisfactorily in the case of $10\text{\,}\mathrm{\char
37\relax}$ noise on the data;

  * •

during the transition from the scalar reconstruction to the vector-valued one,
reconstruction errors can be propagated.

Summarizing, we can say that even though the scalar spline method yields very
good results during the synthetic test case, the vector method should be
preferred if one is interested in a reconstruction of the vector-valued
current. Due to the vector spline method, quality of the reconstruction can be
gained which manifests in the relative NRMSE as well in the deviation plots.
Especially with increasing noise level, the vector spline method is more
robust than the scalar spline method for solving the vector-valued inverse
problem.

Due to very good numerical results throughout our synthetic test cases, we
test the numerical method with real data. To the knowledge of the author, it
is the first time that a reproducing kernel based spline method is applied to
real magneto-electroencephalography data, even though it yielded promising
results for the scalar spline approach in [20].