# Multilinear Kernel Regression and Imputation
via Manifold Learning

Duc Thien Nguyen and Konstantinos Slavakis1 1D. T. Nguyen and K. Slavakis are
with the Department of Information and Communications Engineering, Tokyo
Institute of Technology, Yokohama, 226-8502 Japan (e-mails: {nguyen.t.au,
slavakis.k.aa}@m.titech.ac.jp).

###### Abstract

This paper introduces a novel nonparametric framework for data imputation,
coined multilinear kernel regression and imputation via the manifold
assumption (MultiL-KRIM). Motivated by manifold learning, MultiL-KRIM models
data features as a point cloud located in or close to a user-unknown smooth
manifold embedded in a reproducing kernel Hilbert space. Unlike typical
manifold-learning routes, which seek low-dimensional patterns via regularizers
based on graph-Laplacian matrices, MultiL-KRIM builds instead on the intuitive
concept of tangent spaces to manifolds and incorporates collaboration among
point-cloud neighbors (regressors) directly into the data-modeling term of the
loss function. Multiple kernel functions are allowed to offer robustness and
rich approximation properties, while multiple matrix factors offer low-rank
modeling, integrate dimensionality reduction, and streamline computations with
no need of training data. Two important application domains showcase the
functionality of MultiL-KRIM: time-varying-graph-signal (TVGS) recovery, and
reconstruction of highly accelerated dynamic-magnetic-resonance-imaging (dMRI)
data. Extensive numerical tests on real and synthetic data demonstrate MultiL-
KRIM’s remarkable speedups over its predecessors, and outperformance over
prevalent “shallow” data-imputation techniques, with a more intuitive and
explainable pipeline than deep-image-prior methods.

###### Index Terms:

Imputation, regression, nonparametric, kernel, manifold, graph, MRI

##  I Introduction

Missing data occur frequently across different fields, such as recommender
systems [1], remote sensing [2], sensor networks [3], and accelerated medical
imaging [4, 5], inflicting bias and deficiency on downstream data-analysis
tasks. Amidst a wide range of strategies for data imputation, imputation-by-
regression is one of the most popular frameworks, where measured data act as
regressors while missing data are imputed by the fitted regression model and
measured data [6].

Imputation of missing data is typically viewed as a matrix/tensor completion
problem, since (dynamic) images, or spatio-temporal data of sensor networks
are usually arranged into matrices/tensors. Popular data-imputation approaches
are compressed sensing, e.g., [7], and low-rank methods which are realized via
a decomposition of the data matrix into low-rank matrix factors [8, 9]. On the
other hand, dictionary learning [10, 11, 12, 13] factorizes the data matrix
into an over-complete (“fat”) matrix and a sparse one. Sequential multi-layer
factorization of the data matrix has been also reported [14]. Non-linear
activation function can be also applied to those layers, resulting into deep
matrix factorizations [15, 16], similarly to deep neural networks. Tensors can
be reshaped (“flattened”) into matrices to enjoy the tools of matrix
completion. However, there exist concepts of multilinear rank inherent to
tensors, such as the CP and Tucker rank, yielding CP and Tucker tensor
decompositions [17]. Recently, tensor networks offer generalized tensor
decompositions [18, 19, 20]. All of the aforementioned methods fall under the
umbrella of parametric modeling, where the characteristic is that model
dimensions, that is, inner dimensions of the matrix/tensor factorizations, do
not scale with the number of the observed data.

Another recently popular path of parametric modeling is deep learning (DeepL),
e.g., [21, 22]. Nevertheless, not only is DeepL usually training-data hungry
and computationally heavy, but concerns were also raised in [23] about
instabilities in medical-image reconstruction. Deep image priors (DIP) [24]
and untrained neural networks offer user-defined priors to alleviate the need
for massive training data, and have been used in signal recovery [25, 26] and
accelerated medical imaging [27, 28].

In contrast to the previous parametric methods, nonparametric regression
operates without any priors and statistical assumptions on the data in an
effort to reduce as much as possible the bias inflicted on data modeling by
the user [29]. The price to be paid for this distribution-free approach is
that the dimensions of the nonparametric-regression model scales with the
number of observed data. Within nonparametric-regression methods, the kernel-
based ones stand out thanks to their highly intuitive and explainable
functional forms to approximate involved non-linear functions [30, 31, 32, 33,
34]. The crux is to use a non-linear feature mapping $\varphi(\cdot)$ to map
data from their input space to a high-, even infinite-dimensional feature
space, where linear operations are used to reflect non-linear operations in
the original input space. If the feature space is chosen to be a reproducing
kernel Hilbert space (RKHS) $\mathscr{H}$, with well-known and versatile
properties in approximation theory [35], inner products in $\mathscr{H}$ can
be calculated via simple functional evaluations in the input space (“kernel
trick”) [36].

Aiming at learning methods for high-dimensional data, manifold learning (ManL)
offers a framework where data/features are considered to be points in or close
to a low-dimensional (non-linear) smooth manifold $\mathscr{M}$ embedded in a
high-dimensional ambient space [37]. ManL affords many degrees of freedom in
data modeling because it adopts a minimal number of assumptions; no specifics
of $\mathscr{M}$ are presumed other than its smoothness. ManL identifies low-
dimensional data patterns to effect dimensionality reduction. More
specifically, “local approaches,” including local linear embedding (LLE) [38]
and Laplacian eigenmaps [39], effect dimensionality reduction by preserving
local proximity among data, while “global approaches,” such as ISOMAP [40],
preserve geodesic distances between all pairs of points. While local
approaches may fail to capture long-range dependencies in the point-cloud,
global ones are usually computationally costly. The standard ManL route for
data imputation is to use a graph Laplacian matrix to model the
interdependencies among the rows/columns of the data matrix [41, 42]. This
approach has been widely used in recommender systems [43, 44, 45] to construct
users/items interaction. Due to the focus on “locality” by the graph Laplacian
matrix, the computed interpolation functions appear to be non-continuous at
sampled points, especially in cases of severe under-sampling [46].

Different from parametric modeling and motivated by ManL, this paper develops
a generalization, abbreviated hereafter MultiLinear Kernel Regression and
Imputation via the Manifold assumption (MultiL-KRIM), of the nonparametric
kernel-based approach of [47, 48]. Features are assumed to form a point-cloud
which lies in or close to a smooth manifold $\mathscr{M}$ embedded in an RKHS
$\mathscr{H}$. Unlike popular decomposition-based methods, such as low-rank-
[49, 50], tensor-based [51, 52] and dictionary learning [10, 11, 12, 13],
which promote a “blind decomposition” of the data matrix/tensor, and unlike
standard ManL routes which are based on Laplacian-matrix-based regularization
of a regression loss [41, 42], MultiL-KRIM explores collaborative-filtering
and tangent-space modeling ideas to extract and at the same time incorporate
latent feature geometry directly into data-matrix decompositions. MultiL-KRIM
needs no training data to operate, belongs to the family of nonparametric
regression to reduce its dependence on the probability distribution of the
data as much as possible [29], and offers an explainable and simple geometric
learning paradigm, unlike the majority of DeepL schemes which are based on
perplexed and cascading non-linear function layers.

Unlike its bilinear predecessors BiLMDM [47] and KRIM [48], where
dimensionality-reduction pre-steps are detached from the regression task,
MultiL-KRIM affords multiple matrix factors to allow dimensionality reduction
be employed directly into the regression task by letting its inverse-problem
automatically identify the “optimal low-dimensional” rendition of a kernel
matrix. In addition, MultiL-KRIM leverages its flexible matrix factors for
higher computational efficiency in its inverse-problem algorithmic solution.

To showcase the functionality of MultiL-KRIM, two imputation tasks,
corresponding to different application domains, are considered:  (i) the
recently popular problem of time-varying graph-signal (TVGS) recovery (Section
IV), and  (ii) signal reconstruction in dynamic magnetic resonance imaging
(dMRI) under severe under-sampling (Section V). This manuscript extends the
short paper [53], which considers the dMRI case.  The extensive numerical
tests in Sections IV and V illustrate that MultiL-KRIM outperforms various
state-of-the-art “shallow” parametric and nonparametric methods.

The rest of the paper is organized as follows. Section II describes rigorously
the missing-data imputation problem and briefly introduces state-of-the-art
techniques. Section III presents the MultiL-KRIM data modeling, its inverse
problem and an algorithmic solution. Finally, Section VI summarizes the paper
and provides hints about future research directions. To abide by the Journal’s
page limits for a first paper submission, several of the discussions are moved
into the accompanying supplementary file.

##  II The Imputation-by-Regression Problem and Prior Art

Data are arranged in an $N$-mode tensor $\mathbfscr{Y}$ of size $I_{1}\times
I_{2}\times\cdots\times I_{N}$, where
$\\{I_{n}\\}_{n=1}^{N}\subset\mathbb{N}_{*}$ ($\mathbb{N}_{*}$ is the set of
all positive integers). The observed data take values from the set of all
complex numbers $\mathbb{C}$, while entries of $\mathbfscr{Y}$ which
correspond to missing data are set to be equal to $+\infty$. In other words,
$\mathbfscr{Y}\in(\mathbb{C}\cup\\{+\infty\\})^{I_{1}\times
I_{2}\times\cdots\times I_{N}}$.

Tensor $\mathbfscr{Y}$ can be unfolded/matricized/flattened along mode $k$
into a matrix $\mathbf{Y}_{(k)}$ [17]. For simplicity, assume that
$\mathbfscr{Y}$ is unfolded along mode $N$, and
$\mathbf{Y}\coloneqq\mathbf{Y}_{(N)}\in(\mathbb{C}\cup\\{+\infty\\})^{I_{0}\times
I_{N}}$, where $I_{0}\coloneqq I_{1}I_{2}\cdots I_{N-1}$. Let the $(i,t)$th
entry of $\mathbf{Y}$ be denoted by $y_{it}\coloneqq[\mathbf{Y}]_{it}$. Denote
also the rows and columns of $\mathbf{Y}$ as
$\mathcal{Z}\coloneqq\\{\mathbf{z}_{i}\\}_{i=1}^{I_{0}}$ and
$\mathcal{Y}\coloneqq\\{\mathbf{y}_{t}\\}_{t=1}^{I_{N}}$, respectively; that
is,
$\mathbf{Y}=[\mathbf{z}_{1},\ldots,\mathbf{z}_{I_{0}}]^{\intercal}=[\mathbf{y}_{1},\ldots,\mathbf{y}_{I_{N}}]$,
where $\intercal$ denotes vector/matrix transposition. For rigorous
discussions, let the index set of observed entries
$\Omega\coloneqq\\{(i,t)\in\\{1,\ldots,I_{0}\\}\times\\{1,\ldots,I_{N}\\}\nonscript\,|\allowbreak\nonscript\,\mathopen{}y_{it}\neq+\infty\\}$,
and define the linear sampling mapping
$\mathscr{S}_{\Omega}\colon(\mathbb{C}\cup\\{+\infty\\})^{I_{0}\times
I_{N}}\to\mathbb{C}^{I_{0}\times
I_{N}}\colon\mathbf{Y}\mapsto\mathscr{S}_{\Omega}(\mathbf{Y})$, which operates
entry-wisely as follows:
$[\mathscr{S}_{\Omega}(\mathbf{Y})]_{it}\coloneqq[\mathbf{Y}]_{it}$, if
$(i,t)\in\Omega$, while $[\mathscr{S}_{\Omega}(\mathbf{Y})]_{it}\coloneqq 0$,
if $(i,t)\notin\Omega$.

To cover a wide range of application domains, it is assumed that data
$\mathbf{Y}$, observed in a input-data domain
$\mathscr{D}_{y}\subset(\mathbb{C}\cup\\{+\infty\\})^{I_{0}\times I_{N}}$,
possess also a representation over another domain
$\mathscr{D}_{x}\subset\mathbb{C}^{I_{0}\times I_{N}}$, and that the “link”
between $\mathscr{D}_{y}$ and $\mathscr{D}_{x}$ is a transform mapping
$\mathscr{T}\colon\mathscr{D}_{x}\to\mathscr{D}_{y}$, which is usually
considered to be invertible. In other words, it is assumed that there exists
an $\mathbf{X}\in\mathscr{D}_{x}\subset\mathbb{C}^{I_{0}\times I_{N}}$ such
that (s.t.) $\mathscr{T}(\mathbf{X})$ agrees with
$\mathbf{Y}\in\mathscr{D}_{y}$ on its observed entries at positions $\Omega$.
For example, in the context of Section V, $\mathscr{D}_{y}$ is the k-space
domain and $\mathscr{D}_{x}$ is the image one, with $\mathscr{T}$ being the
two-dimensional discrete Fourier transform (DFT) $\mathscr{F}$ [5]. In cases
where there is no need for the extra domain $\mathscr{D}_{x}$ (see Section
IV), then $\mathscr{D}_{x}\coloneqq\mathscr{D}_{y}$, and $\mathscr{T}$ is set
to be the identity mapping $\operatorname{Id}$.

Typically, imputation-by-regression is achieved by solving the inverse problem

$\displaystyle\min\nolimits_{\mathbf{X}\in\mathbb{C}^{I_{0}\times I_{N}}}{}$
$\displaystyle{}\mathcal{L}(\mathbf{X})+\mathcal{R}(\mathbf{X}){}$ (1a) s.to
$\displaystyle{}\mathscr{S}_{\Omega}(\mathbf{Y})=\mathscr{S}_{\Omega}\mathscr{T}(\mathbf{X})$
(1b) and other constraints,

where $\mathcal{L}(\cdot)$ is a loss that incorporates the user-defined data
model into the design, and $\mathcal{R}(\cdot)$ is the regularizer used to
impose structural priors to the data model. Moreover, (1b) is the “data-
consistency” constraint. A popular strategy for low-rank models is to have
data consistency appear as the “soft” loss
$\mathcal{L}(\mathbf{X})\coloneqq\lVert\mathscr{S}_{\Omega}(\mathbf{Y})-\mathscr{S}_{\Omega}\mathscr{T}(\mathbf{X})\rVert^{2}_{\textnormal{F}}$
in (1a), where $\lVert\cdot\rVert_{\textnormal{F}}$ stands for the Frobenius
norm of a matrix, together with the nuclear-norm regularizer
$\mathcal{R}(\mathbf{X})\coloneqq\lVert\mathbf{X}\rVert_{*}$ [7], instead of
using the “hard” constraint in (1b). Although such a path can be also
addressed by the proposed MultiL-KRIM, the present manuscript keeps the “hard”
constraint in (1b) to allow for $\mathcal{L}(\cdot)$ obtain simpler,
“cleaner,” and more familiar expressions than those after involving the
operator $\mathscr{S}_{\Omega}$ in $\mathcal{L}(\cdot)$.

A popular low-rank data model is expressed by the loss
$\mathcal{L}(\mathbf{X},\mathbf{U},\mathbf{V})=\lVert\mathbf{X}-\mathbf{U}\mathbf{V}\rVert^{2}_{\textnormal{F}}$,
where $\mathbf{U}$ and $\mathbf{V}$ are added as auxiliary variables to the
inverse problem (1), with $\mathbf{U}$ having more rows than columns [8, 9].
In contrast, if $\mathbf{U}$ is over-complete, i.e., more columns than rows,
and $\mathbf{V}$ is sparse, the previous loss yields dictionary learning [11].
More generally, the loss
$\lVert\mathbf{X}-\mathbf{U}_{1}\mathbf{U}_{2}\cdots\mathbf{U}_{Q}\mathbf{V}\rVert_{\textnormal{F}}^{2}$,
where all of $\mathbf{X}$, $\\{\mathbf{U}_{q}\\}_{q=1}^{Q}$ and $\mathbf{V}$
are variables to be optimized, leads to multi-layer matrix factorization (MMF)
[14], which serves as a multilinear representation of the data. With regards
to multilinear algebra, inverse problems can also be formulated directly on
the tensor $\mathbfscr{Y}$, where tensor-decomposition frameworks such as
Tucker and CP decomposition [17], or the more recent tensor networks [18, 19,
20], replace the standard matrix ones. All of the aforementioned parametric
factorization schemes are “blind,” in a sense that all matrix/tensor factors
are variables to be computed, and no data/feature geometry is incorporated in
the data model. Typically, any prior knowledge and additional model
assumptions are realized via the regularizer $\mathcal{R}(\cdot)$.

The nonparametric basis pursuit (NBP) [30] introduces the following data-
modeling loss

$\displaystyle\mathcal{L}(\mathbf{X},\mathbf{B},\mathbf{C})=\lVert\mathbf{X}-\mathbf{K}_{\mathcal{Z}}\mathbf{B}\mathbf{C}\mathbf{K}_{\mathcal{Y}}\rVert^{2}_{\textnormal{F}}\,,$
(2)

where
$[\mathbf{K}_{\mathcal{Z}}]_{i^{\prime}i^{\prime\prime}}=\kappa_{\mathcal{Z}}(\mathbf{z}_{i^{\prime}},\mathbf{z}_{i^{\prime\prime}})$
and
$[\mathbf{K}_{\mathcal{Y}}]_{t^{\prime}t^{\prime\prime}}=\kappa_{\mathcal{Y}}(\mathbf{y}_{t^{\prime}},\mathbf{y}_{t^{\prime\prime}})$
are kernel matrices, generated by the rows $\mathcal{Z}$ and columns
$\mathcal{Y}$, with reproducing kernel functions
$\kappa_{\mathcal{Z}}(\cdot,\cdot)$ and $\kappa_{\mathcal{Y}}(\cdot,\cdot)$
corresponding to some user-defined RKHSs $\mathcal{H}_{\mathcal{Z}}$ and
$\mathcal{H}_{\mathcal{Y}}$, respectively. A very short primer on RKHSs can be
found in Appendix A (supplementary file). The motivation behind (2) is model
$x_{it}\approx\sum_{k=1}^{d}f_{k}(\mathbf{z}_{i})g_{k}(\mathbf{y}_{t})$, where
$f_{k}\coloneqq\sum_{i^{\prime}=1}^{I_{0}}b_{i^{\prime}k}\varphi_{\mathcal{Z}}(\mathbf{z}_{i^{\prime}})\in\mathcal{H}_{\mathcal{Z}}$,
$g_{k}\coloneqq\sum_{t^{\prime}=1}^{I_{N}}c_{kt^{\prime}}\varphi_{\mathcal{Y}}(\mathbf{y}_{t^{\prime}})\in\mathcal{H}_{\mathcal{Y}}$,
with $\varphi_{\mathcal{Z}}$ and $\varphi_{\mathcal{Y}}$ being feature maps
(see Appendix A), and $d\leq\min\\{I_{0},I_{N}\\}$ is a user-defined upper
bound on $\mathbf{X}$’s rank. Matrices $\mathbf{B}$ and $\mathbf{C}$ in (2)
gather all coefficients $\\{b_{i^{\prime}k}\\}$ and $\\{c_{kt^{\prime}}\\}$,
respectively. Similarly to NBP [30], kernel graph learning (KGL) [54] assumes
the loss
$\mathcal{L}(\mathbf{X},\mathbf{G})=\lVert\mathbf{X}-\mathbf{K}_{\mathcal{Z}}\mathbf{G}\mathbf{K}_{\mathcal{Y}}\rVert^{2}_{\textnormal{F}}$,
which yields (2) by decomposing $\mathbf{G}$ as the product $\mathbf{BC}$ of
low-rank matrices. Moreover, kernel regression over graphs (KRG) [55] adopts
the loss
$\mathcal{L}(\mathbf{X},\mathbf{H})=\lVert\mathbf{X}-\mathbf{H}\mathbf{K}_{\mathcal{Y}}\rVert^{2}_{\textnormal{F}}$,
which yields NBP and KGL by decomposing $\mathbf{H}$ into
$\mathbf{K}_{\mathcal{Z}}\mathbf{BC}$ and
$\mathbf{K}_{\mathcal{Z}}\mathbf{G}$, respectively. To avoid alternating
minimization for solving for the bilinear (2), [33] proposed the loss
$\mathcal{L}(\mathbf{X},\bm{\gamma})=\lVert\operatorname{vec}(\mathbf{X})-\mathbf{K}_{\otimes}\bm{\gamma}\rVert^{2}_{\textnormal{F}}$,
where $\mathbf{K}_{\otimes}$ is the Kronecker product
$\mathbf{K}_{\mathcal{Y}}\otimes\mathbf{K}_{\mathcal{Z}}$ and
$\operatorname{vec}(\mathbf{X})$ is the vectorization of $\mathbf{X}$.
Although [33] offers an efficient closed-form solution, its data modeling
lacks low-rankness and bilinearity. Along the lines of KRG [55] and motivated
by LLE [38], [56] introduces
$\mathcal{L}(\mathbf{X},\mathbf{H},\mathbf{F},\mathbf{Q})=\lVert\mathbf{X}-\mathbf{H}\mathbf{K}_{F}^{2}\rVert_{\textnormal{F}}+\lVert\mathbf{F}-\mathbf{Q}\mathbf{K}_{\mathcal{Y}}\rVert^{2}_{\textnormal{F}}$,
where the kernel matrix $\mathbf{K}_{F}$ is defined by the low-dimensional
representation $\mathbf{F}$ of $\mathbf{X}$. This method requires two
optimization steps per iteration, first on $(\mathbf{H},\mathbf{Q})$, and then
on $(\mathbf{X},\mathbf{F})$.

Notice that kernel matrices
$\mathbf{K}_{\mathcal{Z}},\mathbf{K}_{\mathcal{Y}}$ are square, with their
dimensions scaling with the dimensions of $\mathbf{X}$. This causes the
scaling of the dimensions of the parameter matrices
$\mathbf{B},\mathbf{C},\mathbf{G}$, and $\mathbf{H}$ with the dimensions of
the data, which is characteristic of nonparametric methods, but at the same
time raises significant computational obstacles in cases of high-dimensional
data. In fact, the CPU setting, used in this study for numerical validations,
was not able to carry through the NBP computations for the high-dimensional
data of Section V. Moreover, all approximating functions $\\{f_{k}\\}$ and
$\\{g_{k}\\}$ in NBP, KGL, and KRG are only required to belong to some RKHSs,
without any consideration of more elaborate geometric and ManL arguments.

##  III Data Modeling by MultiL-KRIM

Step 1 Identify navigator data $\check{\mathbf{Y}}_{\textnormal{nav}}$ Step 2
Select landmark points $\check{\mathbf{L}}$ Step 3 Compute kernel matrix
$\mathbfcal{K}$ in (17) (supplementary file) Step 4 Design and solve the
inverse problem (11) Figure 1: Pipeline of MultiL-KRIM
$\varphi(\mathbfit{l}_{k_{1}})$$\varphi(\mathbfit{l}_{k_{2}})$$\varphi(\mathbfit{l}_{k_{3}})$$\varphi(\check{\bm{\mu}}_{t})$$\mathscr{M}$$\mathscr{H}$
Figure 2: The manifold-learning (ManL) modeling assumption: points
$\\{\varphi(\mathbfit{l}_{k})\\}_{k=1}^{N_{\mathit{l}}}$ lie into or close to
an unknown-to-the-user smooth manifold $\mathscr{M}$ which is embedded into an
ambient RKHS $\mathscr{H}$. The “collaborative-filtering” modeling assumption:
only a few points, $\\{\varphi(\mathbfit{l}_{k_{n}})\\}_{n=1}^{3}$ here,
collaborate affinely to approximate $\varphi(\check{\bm{\mu}}_{t})$. The small
number of “collaborating neighbors” implies a low-dimensional structure. All
affine combinations of $\\{\varphi(\mathbfit{l}_{k_{n}})\\}_{n=1}^{3}$ define
the approximating “linear patch” (affine hull, gray-colored plane), which
mimics the concept of a tangent space to $\mathscr{M}$.

As a first step (see Figure 1), MultiL-KRIM extracts geometric information
from a subset of the observed data, coined “navigator data” and denoted
hereafter as
$\check{\mathbf{Y}}_{\textnormal{nav}}\coloneqq[\check{\mathbf{y}}_{1}^{\textnormal{nav}},\ldots,\check{\mathbf{y}}_{N_{\text{nav}}}^{\textnormal{nav}}]\in\mathbb{C}^{\nu\times
N_{\text{nav}}}$, for some $\nu,N_{\text{nav}}\in\mathbb{N}_{*}$. Their
specific definition depends on the application domain. For example, navigator
data can be vectorized patches from the heavily sampled region of an image,
or, in a network with time-series data, they can be samples of time-series
collected from a subset of the network nodes. This approach of extracting
geometric information from the data themselves, without imposing any user-
defined statistical priors, and incorporating the extracted information into
MultiL-KRIM’s inverse problem differs from “blind” decomposition methods,
e.g., [49, 50, 51, 52, 10, 11, 12, 13] (see also Section I), where no
geometric information from the data is infused into their inverse problems.

As the number of navigator data points $N_{\text{nav}}$ grows for large
datasets, to avoid computational complexity issues, a subset
$\\{\mathbfit{l}_{k}\\}_{k=1}^{N_{\mathit{l}}}$, coined
landmark/representative points, with $N_{\mathit{l}}\leq N_{\text{nav}}$, is
selected from
$\\{\check{\mathbf{y}}^{\textnormal{nav}}_{k}\\}_{k=1}^{N_{\text{nav}}}$ (see
Figure 1). To avoid interrupting the flow of the discussion, several
strategies to form navigator data $\check{\mathbf{Y}}_{\textnormal{nav}}$ and
to select landmark points $\\{\mathbfit{l}_{k}\\}_{k=1}^{N_{\mathit{l}}}$ are
deferred to Section IV-C. Let
$\check{\mathbf{L}}\coloneqq[\mathbfit{l}_{1},\mathbfit{l}_{2},\ldots,\mathbfit{l}_{N_{\mathit{l}}}]\in\mathbb{C}^{\nu\times
N_{\mathit{l}}}$. Motivated by ManL, this study assumes that there is some
latent geometry hidden within
$\\{\check{\mathbf{y}}^{\textnormal{nav}}_{k}\\}_{k=1}^{N_{\text{nav}}}$ and
consequently $\\{\mathbfit{l}_{k}\\}_{k=1}^{N_{\mathit{l}}}$, which manifests
itself in an appropriately user-defined feature space $\mathscr{H}$. Here,
$\mathscr{H}$ is taken to be an RKHS, associated with a reproducing kernel
$\kappa(\cdot,\cdot)\colon\mathbb{C}^{\nu}\times\mathbb{C}^{\nu}\to\mathbb{R}$,
and a feature mapping $\varphi\colon\mathbb{C}^{\nu}\to\mathscr{H}$ (see
Appendix A in the supplementary file). By a slight abuse of notation, let
$\bm{\Phi}(\check{\mathbf{L}})\coloneqq[\varphi(\mathbfit{l}_{1}),\ldots,\varphi(\mathbfit{l}_{N_{\mathit{l}}})]$,
where $\varphi(\mathbfit{l}_{k})\in\mathscr{H}$, $\forall
k\in\\{1,\ldots,N_{\mathit{l}}\\}$ (see Figure 2).

The $(i,t)$th entry $x_{it}$ of $\mathbf{X}$ is modeled as

$\displaystyle[\mathbf{X}]_{it}=x_{it}$ $\displaystyle{}\approx{}$
$\displaystyle f_{i}(\check{\bm{\mu}}_{t})=\langle
f_{i}\mathop{}|\mathop{}\varphi(\check{\bm{\mu}}_{t})\rangle_{\mathscr{H}}\,,$
(3)

where $f_{i}(\cdot):\mathbb{C}^{\nu}\to\mathbb{C}$ is an unknown non-linear
function in the RKHS functional space $\mathscr{H}$, $\check{\bm{\mu}}_{t}$ is
an unknown vector in $\mathbb{C}^{\nu}$, and the latter part of (3) is because
of the reproducing property of $\mathscr{H}$ [35]. It is assumed that $f_{i}$
and $\varphi(\check{\bm{\mu}}_{t})$ belong to the linear span of
$\\{\varphi(\mathbfit{l}_{k})\\}_{k=1}^{N_{\mathit{l}}}$, i.e., there exist
complex-valued $N_{\mathit{l}}\times 1$ vectors
$\mathbf{d}_{i}\coloneqq[d_{i1},\ldots,d_{iN_{\mathit{l}}}]^{\intercal}$ and
$\mathbf{b}_{t}\coloneqq[b_{1t},\ldots,b_{N_{\mathit{l}}t}]^{\intercal}$ s.t.

$\displaystyle f_{i}$
$\displaystyle=\sum\nolimits_{k=1}^{N_{\mathit{l}}}d_{ik}\varphi(\mathbfit{l}_{k})=\bm{\Phi}(\check{\mathbf{L}})\mathbf{d}_{i}\,,$
(4a) $\displaystyle\varphi(\check{\bm{\mu}}_{t})$
$\displaystyle=\sum\nolimits_{k=1}^{N_{\mathit{l}}}b_{kt}\varphi(\mathbfit{l}_{k})=\bm{\Phi}(\check{\mathbf{L}})\mathbf{b}_{t}\,.$
(4b)

Altogether, by back-substitutions into (3),

$\displaystyle
x_{it}\approx\langle\bm{\Phi}(\check{\mathbf{L}})\mathbf{d}_{i}\mathop{}|\mathop{}\bm{\Phi}(\check{\mathbf{L}})\mathbf{b}_{t}\rangle_{\mathscr{H}}=\mathbf{d}_{i}^{\mathsf{H}}\mathbf{K}\mathbf{b}_{t}\,,$
(5)

where $\mathbf{K}$ is the complex-valued $N_{\mathit{l}}\times N_{\mathit{l}}$
matrix whose $(k,k^{\prime})$-th entry is equal to
$\langle\varphi(\mathbfit{l}_{k})\mathop{}|\mathop{}\varphi(\mathbfit{l}_{k^{\prime}})\rangle_{\mathscr{H}}=\kappa(\mathbfit{l}_{k},\mathbfit{l}_{k^{\prime}})$,
and superscript $\mathsf{H}$ stands for complex conjugate transposition. To
offer a simple data model and compact notations, if
$\mathbf{D}\coloneqq[\mathbf{d}_{1},\ldots,\mathbf{d}_{I_{0}}]^{{\mathsf{H}}}\in\mathbb{C}^{I_{0}\times
N_{\mathit{l}}}$ and
$\mathbf{B}\coloneqq[\mathbf{b}_{1},\ldots,\mathbf{b}_{I_{N}}]\in\mathbb{C}^{N_{\mathit{l}}\times
I_{N}}$ are the parameters that need to be identified, then data are modeled
as the bilinear form

$\displaystyle\mathbf{X}\approx\mathbf{D}\mathbf{K}\mathbf{B}\,.$ (6)

So far, (4) follow the standard subspace-modeling assumption of kernel
methods, met also in the prior nonparametric designs NBP [30], KRG [55] and
KGL [54] (see also the discussion after (2)): both $f_{i}$ and
$\varphi(\check{\bm{\mu}}_{t})$ are assumed to belong to a linear subspace
spanned by $\\{\varphi(\mathbfit{l}_{k})\\}_{k=1}^{N_{\mathit{l}}}$. However,
inspired by ManL and [47, 48], and in quest for low-dimensional structures,
MultiL-KRIM takes the standard kernel modeling assumptions a few steps
further: points $\\{\varphi(\mathbfit{l}_{k})\\}_{k=1}^{N_{\mathit{l}}}$ are
assumed to lie into or close to an unknown-to-the-user smooth manifold
$\mathscr{M}$ which is embedded into an ambient RKHS $\mathscr{H}$ (see Figure
2). Furthermore, a “collaborative-filtering” modeling flavor is added to the
design: only a few points within
$\\{\varphi(\mathbfit{l}_{k})\\}_{k=1}^{N_{\mathit{l}}}$
($\\{\varphi(\mathbfit{l}_{k_{n}})\\}_{n=1}^{3}$ in Figure 2) collaborate
affinely to approximate $\varphi(\check{\bm{\mu}}_{t})$. Affine combinations,
expressed via $\mathbf{1}_{N_{\mathit{l}}}^{{\mathsf{H}}}\mathbf{b}_{t}=1$,
where $\mathbf{1}_{N_{\mathit{l}}}$ is the $N_{\mathit{l}}\times 1$ all-one
vector, are motivated by the fundamental concept of tangent spaces (“linear
patches”, simply put) to smooth manifolds [57]. Moreover, the assumption that
only a small number of “collaborating neighbors” is sufficient to approximate
$\varphi(\check{\bm{\mu}}_{t})$, expressed explicitly by considering a sparse
vector $\mathbf{b}_{t}$, is motivated by the desire to identify low-
dimensional structures within
$\\{\varphi(\mathbfit{l}_{k})\\}_{k=1}^{N_{\mathit{l}}}$ and by the success of
sparse approximations in low-dimensional data modeling.

Choosing $\kappa(\cdot,\cdot)$ to define $\mathbf{K}$ requires the intensive
tasks of cross validation and fine tuning. A popular and straightforward way
to overcome such tasks and to generalize (6), implemented also in [48], is via
multiple kernels:

$\displaystyle\mathbf{X}\approx\sum\nolimits_{m=1}^{M}\mathbf{D}_{m}\mathbf{K}_{m}\mathbf{B}_{m}\,,$
(7)

with user-defined reproducing kernels
$\\{\kappa_{m}(\cdot,\cdot)\\}_{m=1}^{M}$, and thus kernel matrices
$\\{\mathbf{K}_{m}\\}_{m=1}^{M}$, complex-valued $I_{0}\times N_{\mathit{l}}$
matrices $\\{\mathbf{D}_{m}\\}_{m=1}^{M}$, and $N_{\mathit{l}}\times I_{N}$
matrices $\\{\mathbf{B}_{m}\\}_{m=1}^{M}$.

The larger the cardinality $N_{\text{nav}}$ of the navigator data, the larger
the number $N_{\mathit{l}}$ of landmark points may become, justifying the
nonparametric character of the present design. Likewise, the larger the
dimensionality $I_{0}$ of the data $\mathbf{Y}$, the larger the number of rows
of $\mathbf{D}_{m}$s in (7). In other words, for large and high-dimensional
datasets (see Section V), the size $I_{0}\times N_{\mathit{l}}$ of
$\mathbf{D}_{m}$s may become unbearably large. To reduce computational
complexity and effect low-rank structure on the approximation, KRIM [48]
compressed each $N_{\mathit{l}}\times N_{\mathit{l}}$ matrix $\mathbf{K}_{m}$
into a low-dimensional $d\times N_{\mathit{l}}$ ($d\ll N_{\mathit{l}}$)
rendition $\check{\mathbf{K}}_{m}$ to yield the model
$\mathbf{X}\approx\sum_{m=1}^{M}\check{\mathbf{D}}_{m}\check{\mathbf{K}}_{m}\mathbf{B}_{m}$,
where $\\{\check{\mathbf{D}}_{m}\\}_{m=1}^{M}$ are now low-rank $I_{0}\times
d$ matrices. Matrix $\check{\mathbf{K}}_{m}$ was computed from
$\mathbf{K}_{m}$ via a dimensionality-reduction module (pre-step) in [48],
inspired by LLE [38]. However, such a pre-step inflicts the following
drawbacks on KRIM’s design:  (i) with large numbers $M$ of kernels and
$N_{\mathit{l}}$ of landmark points, the LLE-inspired pre-step raises
computational obstacles, while fine-tuning its hyperparameters consumes extra
effort; and  (ii) if not carefully fine-tuned, the error from compressing
$\mathbf{K}_{m}$ into $\check{\mathbf{K}}_{m}$ may propagate to the next stage
of KRIM’s pipeline.

The present MultiL-KRIM avoids the previous KRIM’s drawbacks by the following
multilinear factorization,

$\displaystyle\mathbf{X}$
$\displaystyle\approx\sum\nolimits_{m=1}^{M}\mathbf{D}_{m}^{(1)}\mathbf{D}_{m}^{(2)}\cdots\mathbf{D}_{m}^{(Q)}\mathbf{K}_{m}\mathbf{B}_{m}\,,$
(8)
$\displaystyle\stackunder{\fcolorbox{black}{tableSingle}{\makebox[6\matfield]{\addstackgap[.5\tmplength]{\box
0}}}}{\scriptstyle(I_{0}\times I_{N})}$
$\displaystyle\approx\sum_{m=1}^{M}{\stackunder{\fcolorbox{black}{tableSingle}{\makebox[3\matfield]{\addstackgap[.5\tmplength]{\box
0}}}}{\scriptstyle(I_{0}\times d_{1})}}\ \cdots\raisebox{19.5pt}{
\stackunder{\fcolorbox{black}{tableSingle}{\makebox[4\matfield]{\addstackgap[.5\tmplength]{\box
0}}}}{\scriptstyle(d_{Q-1}\times N_{\mathit{l}})}}\raisebox{15.0pt}{
\stackunder{\fcolorbox{black}{white}{\makebox[4\matfield]{\addstackgap[.5\tmplength]{\box
0}}}}{\scriptstyle(N_{\mathit{l}}\times N_{\mathit{l}})}}\raisebox{15.0pt}{
\stackunder{\fcolorbox{black}{tableSingle}{\makebox[6\matfield]{\addstackgap[.5\tmplength]{\box
0}}}}{\scriptstyle(N_{\mathit{l}}\times I_{N})}}$

where $\mathbf{D}_{m}^{(q)}\in\mathbb{C}^{d_{q-1}\times d_{q}}$,
$d_{0}\coloneqq I_{0}$, $d_{Q}\coloneqq N_{\mathit{l}}$, and the inner matrix
dimensions $\\{d_{q}\\}_{q=1}^{Q-1}$ are user defined. Notice that for $Q=2$,
the term $\mathbf{D}_{m}^{(2)}\mathbf{K}_{m}$ may be seen as the low-
dimensional $\check{\mathbf{K}}_{m}$ in [48]. Nonetheless,
$\\{\mathbf{D}_{m}^{(q)}\\}_{(q,m)}$ are identified during a single-stage
learning task, avoiding any KRIM-like pre-steps with their hyperparameter
tuning and errors. It is worth stressing here that in the special case where
$M=1$, $\mathbf{K}_{1}\coloneqq\mathbf{I}_{N_{\mathit{l}}}$, with
$\mathbf{I}_{N_{\mathit{l}}}$ denoting the identity matrix, and with no affine
and sparsity constraints on $\mathbf{B}_{1}$, (8) yields the latent-geometry-
agnostic MMF [14]. See Section IV-D for numerical tests of MultiL-KRIM vs. MMF
on real-world data.

It can be verified that $Q>1$ in (8) may offer substantial computational
savings compared to the $Q=1$ case. Indeed, the number of unknowns to be
identified in (8) for $Q>1$ is

$\displaystyle
N_{Q>1}=M\left(\sum\nolimits_{q=1}^{Q}d_{q-1}d_{q}+I_{N}N_{\mathit{l}}\right)\,,$
(9)

as opposed to

$\displaystyle N_{Q=1}=M(I_{0}N_{\mathit{l}}+I_{N}N_{\mathit{l}})$ (10)

in the $Q=1$ case. If $I_{0}$ and $N_{\mathit{l}}$ are large, then
$\\{d_{q}\\}_{q=1}^{Q-1}$ can be chosen so that $N_{Q>1}\ll N_{Q=1}$ (see
Table IX (supplementary file)).

By constructing the “supermatrices”
$\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q}\,,\mathbfcal{K}$, and $\mathbfcal{B}$,
which are described in Appendix B, the previous considerations yield the
compact data model (16) (supplementary file). Under the previous definitions,
the generic MultiL-KRIM inverse problem becomes

$\displaystyle\min_{(\mathbf{X},\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q},\mathbfcal{B})}{}\
$
$\displaystyle{}\tfrac{1}{2}\lVert\mathbf{X}-\mathbfcal{D}_{1}\mathbfcal{D}_{2}\cdots\mathbfcal{D}_{Q}\mathbfcal{K}\mathbfcal{B}\rVert^{2}_{\textnormal{F}}$
$\displaystyle+\underbrace{\lambda_{1}\lVert\mathbfcal{B}\rVert_{1}+\mathcal{R}_{2}(\mathbf{X},\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q})}_{\mathcal{R}(\mathbf{X},\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q},\mathbfcal{B})}$
(11a) s.to
$\displaystyle{}\mathscr{S}_{\Omega}(\mathbf{Y})=\mathscr{S}_{\Omega}\mathscr{T}(\mathbf{X})\,,$
(11b)
$\displaystyle{}\mathbf{1}_{N_{\mathit{l}}}^{{\mathsf{H}}}\mathbf{B}_{m}=\mathbf{1}_{I_{N}}^{{\mathsf{H}}}\,,\forall
m\in\\{1,\ldots,M\\}\,,$ (11c) and other constraints, (11d)

where the $\ell_{1}$-norm
$\mathcal{R}_{1}(\cdot)\coloneqq\lambda_{1}\lVert\cdot\rVert_{1}$, used here
to impose sparsity on $\mathbfcal{B}$, and $\mathcal{R}_{2}(\cdot)$, used to
enforce prior information on the model, constitute the generic regularizer
$\mathcal{R}\coloneqq\mathcal{R}_{1}+\mathcal{R}_{2}$ of (1a). Constraint
(11b) states the requirement for data consistency over the index set $\Omega$
of observed data values, (11c) stands for the aforementioned ManL-inspired
affine combinations, and (11d) is used to conveniently include any potential
application-domain specific constraint.

The pipeline of MultiL-KRIM framework is illustrated in Figure 1. To validate
the proposed framework, the following sections will demonstrate two important
applications, will provide solutions to specific instances of the generic
inverse problem (11), and will leverage extensive numerical tests to take a
thorough look at the pipeline of Figure 1.

##  IV Application: Time-Varying Graph-Signal Recovery

This section begins with an overview of the problem of time-varying graph-
signal (TVGS) recovery, followed by a short review of several state-of-the-art
prior works in Section IV-A. Section IV-B formulates and provides solutions to
a special instance of the MultiL-KRIM’s inverse problem (11). Finally, Section
IV-C compares empirically MultiL-KRIM with several state-of-the-art methods,
and investigates several variations of MultiL-KRIM.

###  IV-A TVGS recovery

Graph signal processing (GSP) is an emerging field at the intersection of
graph theory and signal processing [58, 59]. A graph is an abstract
mathematical object equipped with nodes/vertices and edges, where nodes
correspond to entities such as sensors of a physical network or data points,
and edges represent the connections/relationships between these entities. In
GSP, each node is annotated with a (time-varying) signal, so that all signals
gathered across all nodes constitute the so-called graph signal.

To be more specific, a graph is denoted by $G=(\mathcal{V},\mathcal{E})$,
where $\mathcal{V}\coloneqq\\{1,2,\ldots,I_{0}\\}$ represents the set of
$I_{0}$ nodes, and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the
set of edges. The topology of the graph is described by the adjacency matrix
$\mathbf{W}\in\mathbb{R}^{I_{0}\times I_{0}}$, whose $(i,j)$th entry
$w_{ij}\in\mathbb{R}_{+}$ denotes the edge weight between nodes $i$ and $j$.
In this study, $G$ is considered to be undirected, thus, $\mathbf{W}$ is
symmetric. The graph Laplacian $I_{0}\times I_{0}$ matrix is defined as
$\mathbf{L}\coloneqq\operatorname{diag}(\mathbf{W}\mathbf{1}_{I_{0}})-\mathbf{W}$,
where $\mathbf{1}_{I_{0}}$ denotes the $I_{0}\times 1$ all-one vector, and
$\operatorname{diag}(\mathbf{W}\mathbf{1}_{I_{0}})$ stands for the
$I_{0}\times I_{0}$ diagonal matrix whose diagonal entries are the
corresponding entries of the $I_{0}\times 1$ vector
$\mathbf{W}\mathbf{1}_{I_{0}}$. The graph signal is denoted by
$\mathbf{Y}=[\mathbf{y}_{1},\mathbf{y}_{2},\ldots,\mathbf{y}_{I_{N}}]\in\mathbb{R}^{I_{0}\times
I_{N}}$, where the columns
$\mathcal{Y}\coloneqq\\{\mathbf{y}_{t}\\}_{t=1}^{I_{N}}$ are the graph-signal
“snapshots” at $I_{N}$ time points, while rows
$\mathcal{Z}\coloneqq\\{\mathbf{z}_{i}\\}_{i=1}^{I_{0}}$ of $\mathbf{Y}$ are
the time profiles of the $I_{0}$ signals which annotate the nodes of the
graph.

The analysis of graph signals is pivotal in fields such as recommender systems
[60], power grids [61], sensor [62], urban [63], and vehicular networks [64].
Due to reasons like users’ privacy, sensors fault, or resources conservation,
such application domains often encounter the problem of missing data.
Naturally, there have been numerous studies to address the problem, most of
which leverage low-rank and smoothness assumptions on $\mathbf{Y}$. Typically,
it is assumed that $\mathbf{Y}$ is low rank due to dependencies among time
points and/or nodes [65]. On the other hand, there are studies which adopt the
graph Laplacian operator $\mathbf{L}$ [58] to impose graph-wide (a.k.a.
spatial) smoothness, i.e., neighboring nodes share similar signal attributes
[66, 67]. For example, the regularization function $\mathcal{R}(\cdot)$ of
both KGL and KRG (see Section II) comprises a classical Tikhonov regularizer
and an $\mathbf{L}$-based quadratic form. Instead, NBP uses only Tikhonov
regularization for $\mathcal{R}(\cdot)$. Although KGL and KRG jointly recover
$\mathbf{Y}$ and the underlying graph topology $\mathbf{W}$ or $\mathbf{L}$,
since this study emphasizes on TVGS recovery, $\mathbf{L}$ is assumed to be
given and fixed in the sequel. Another kernel-based regression method can be
also found in [68], where kernel ridge regression is used to develop a kernel
Kalman-filtering (KKF) algorithm for TVGS recovery, with an $\mathbf{L}$-based
regularizer $\mathcal{R}(\cdot)$. Furthermore, [69, 70] introduce differential
smoothness and impose low-rank structure by nuclear-norm regularizers. Low-
rank-cognizant matrix factorizations have been also studied extensively [71,
72, 73, 74, 75].

By using the one-step difference operator $\mathbf{\Delta}$ [71], study [69]
introduces the spatio-temporal-smoothness regularizer
$\operatorname{tr}(\mathbf{\Delta}^{\intercal}\mathbf{X}^{\intercal}\mathbf{L}\mathbf{X}\mathbf{\Delta})$
in the place of $\mathcal{R}(\cdot)$ in (1), where $\operatorname{tr}(\cdot)$
stands for the trace operator. Specifically, $\mathbf{\Delta}$ is an
$I_{N}\times(I_{N}-1)$ matrix, s.t.
$[\mathbf{\Delta}]_{i,i}=-1,[\mathbf{\Delta}]_{i+1,i}=1$ for all
$i=1,\ldots,I_{N}-1$, while other entries are zeros. Subsequently, [70]
considered an additional low-rank constraint in the form of a nuclear-norm
regularizer, which was shown to outperform other popular regularizers such as
spatial $\operatorname{tr}(\mathbf{X}^{\intercal}\mathbf{LX})$ and temporal
$\lVert\mathbf{X}\mathbf{\Delta}\rVert_{\textnormal{F}}^{2}$ smoothness [76].
On the other hand, [77] introduces the following Sobolev-smoothness
regularizer
$\operatorname{tr}(\mathbf{\Delta}^{\intercal}\mathbf{X}^{\intercal}(\mathbf{L}+\epsilon\mathbf{I}_{I_{0}})^{\beta}\mathbf{X}\mathbf{\Delta})$,
where $\epsilon,\beta\in\mathbb{R}_{++}$ are user-defined parameters. For
convenience, denote
$\mathbf{L}_{\epsilon}^{\beta}\coloneqq(\mathbf{L}+\epsilon\mathbf{I}_{I_{0}})^{\beta}$.
Works [69, 70, 77] use a “soft” data-consistency loss, i.e., the quadratic
loss
$(1/2)\lVert\mathscr{S}_{\Omega}(\mathbf{X})-\mathscr{S}_{\Omega}(\mathbf{Y})\rVert_{\textnormal{F}}^{2}$
is used instead of the “hard” constraint
$\mathscr{S}_{\Omega}(\mathbf{X})=\mathscr{S}_{\Omega}(\mathbf{Y})$.

The extensive numerical tests of the present study have shown that the
Sobolev-smoothness regularizer is currently offering state-of-the-art
performance in TVGS recovery. Consequently, to ensure fairness during
numerical tests, all competing methods, as well as the proposed MultiL-KRIM,
will be equipped hereafter with the Sobolev-smoothness regularizer.

###  IV-B The inverse problem

Algorithm 1 Solving MultiL-KRIM’s inverse problem

0: Limit point $\hat{\mathbf{X}}^{(*)}$ of sequence
$(\hat{\mathbf{X}}^{(n)})_{n\in\mathbb{N}}$.

1: Fix $\hat{\mathbfcal{O}}^{(0)}$, $\gamma_{0}\in(0,1]$, and $\zeta\in(0,1)$.

2: while $n\geq 0$ do

3: Available are $\hat{\mathbfcal{O}}^{(n)}$ ((13) for TVGS or (20) for dMRI)
and $\gamma_{n}$.

4: $\gamma_{n+1}\coloneqq\gamma_{n}(1-\zeta\gamma_{n})$.

5: Solve sub-tasks (14) for TVGS or (21) for dMRI.

6:
$\hat{\mathbfcal{O}}^{(n+1)}\coloneqq\gamma_{n+1}\hat{\mathbfcal{O}}^{(n+1/2)}+(1-\gamma_{n+1})\hat{\mathbfcal{O}}^{(n)}$.

7: Set $n\leftarrow n+1$ and go to step 2.

8: end while

The inverse problem for TVGS recovery becomes a special instance of the
generic (11) one. More specifically,

$\displaystyle\min_{(\mathbf{X},\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q},\mathbfcal{B})}{}\
$
$\displaystyle{}\tfrac{1}{2}\lVert\mathbf{X}-\mathbfcal{D}_{1}\mathbfcal{D}_{2}\cdots\mathbfcal{D}_{Q}\mathbfcal{K}\mathbfcal{B}\rVert^{2}_{\textnormal{F}}+\underbrace{\lambda_{1}\lVert\mathbfcal{B}\rVert_{1}}_{\mathcal{R}_{1}(\mathbfcal{B})}$
$\displaystyle+\underbrace{\tfrac{\lambda_{L}}{2}\operatorname{tr}(\mathbf{\Delta}^{\intercal}\mathbf{X}^{\intercal}\mathbf{L}_{\epsilon}^{\beta}\mathbf{X}\mathbf{\Delta})+\tfrac{\lambda_{2}}{2}\sum_{q=1}^{Q}\lVert\mathbfcal{D}_{q}\rVert_{\textnormal{F}}^{2}}_{\mathcal{R}_{2}(\mathbf{X},\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q})}$
(12a) s.to
$\displaystyle\mathscr{S}_{\Omega}(\mathbf{Y})=\mathscr{S}_{\Omega}(\mathbf{X})\,,$
(12b)
$\displaystyle\mathbf{1}_{N_{\mathit{l}}}^{{\mathsf{H}}}\mathbf{B}_{m}=\mathbf{1}_{I_{N}}^{{\mathsf{H}}}\,,\forall
m\in\\{1,\ldots,M\\}\,,$ (12c)

where the generic regularizer in (11) becomes here
$\mathcal{R}(\mathbf{X},\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q},\mathbfcal{B})\coloneqq\mathcal{R}_{1}(\mathbfcal{B})+\mathcal{R}_{2}(\mathbf{X},\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q})$.
In the current TVGS setting, the transform mapping $\mathscr{T}$ in (11b)
takes the form of the identity operator $\operatorname{Id}$. The loss function
is non-convex, and to guarantee convergence to a critical point, the
successive-convex-approximation framework of [78] is utilized. The algorithm
is summarized in Algorithm 1, where the following tuple of estimates $\forall
n\in\mathbb{N},\forall k\in\\{0,1\\}$,

$\displaystyle\hat{\mathbfcal{O}}^{(n+k/2)}\coloneqq({}$
$\displaystyle{}\hat{\mathbf{X}}^{(n+k/2)},\hat{\mathbfcal{D}}_{1}^{(n+k/2)},$
$\displaystyle\ldots,\hat{\mathbfcal{D}}_{Q}^{(n+k/2)},\hat{\mathbfcal{B}}^{(n+k/2)})\,,$
(13)

is recursively updated via the following convex sub-tasks which need to be
solved at every iteration $n$:

$\displaystyle\hat{\mathbf{X}}^{(n+1/2)}$
$\displaystyle\coloneqq\arg\min_{{\mathbf{X}}}{}$ $\displaystyle{}\
\tfrac{1}{2}\lVert{\mathbf{X}}-\hat{\mathbfcal{D}}_{1}^{(n)}\cdots\hat{\mathbfcal{D}}_{Q}^{(n)}\mathbfcal{K}\hat{\mathbfcal{B}}^{(n)}\rVert_{\textnormal{F}}^{2}$
$\displaystyle{}\
+\tfrac{\lambda_{L}}{2}\operatorname{tr}(\mathbf{\Delta}^{\intercal}\mathbf{X}^{\intercal}\mathbf{L}_{\epsilon}^{\beta}\mathbf{X}\mathbf{\Delta})$
$\displaystyle{}\
+\tfrac{\tau_{X}}{2}\lVert\mathbf{X}-\hat{\mathbf{X}}^{(n)}\rVert_{\textnormal{F}}^{2}$
s.to $\displaystyle{}\
\mathscr{S}_{\Omega}(\mathbf{Y})=\mathscr{S}_{\Omega}({\mathbf{X}})\,,$ (14a)
$\displaystyle\hat{\mathbfcal{D}}_{q}^{(n+1/2)}{}$
$\displaystyle{}\coloneqq\arg\min_{\mathbfcal{D}_{q}}{}$ $\displaystyle{}\
\tfrac{1}{2}\lVert\hat{\mathbf{X}}^{(n)}-\hat{\mathbfcal{D}}_{1}^{(n)}\cdots\mathbfcal{D}_{q}\cdots\hat{\mathbfcal{D}}_{Q}^{(n)}\mathbfcal{K}\hat{\mathbfcal{B}}^{(n)}\rVert_{\textnormal{F}}^{2}$
$\displaystyle{}\
+\tfrac{\lambda_{2}}{2}\lVert\mathbfcal{D}_{q}\rVert_{\textnormal{F}}^{2}+\tfrac{\tau_{D}}{2}\lVert\mathbfcal{D}_{q}-\hat{\mathbfcal{D}}_{q}^{(n)}\rVert_{\textnormal{F}}^{2}$
s.to $\displaystyle{}\ \mathbfcal{D}_{q}\ \textnormal{is block diagonal}\
\forall q\in\\{2,\ldots,Q\\}\,,$ (14b)
$\displaystyle\hat{\mathbfcal{B}}^{(n+1/2)}$
$\displaystyle\coloneqq\arg\min_{\mathbfcal{B}}{}$ $\displaystyle{}\
\tfrac{1}{2}\lVert\hat{\mathbf{X}}^{(n)}-\hat{\mathbfcal{D}}_{1}^{(n)}\cdots\hat{\mathbfcal{D}}_{Q}^{(n)}\mathbfcal{K}\mathbfcal{B}\rVert_{\textnormal{F}}^{2}$
$\displaystyle{}\
+\lambda_{1}\lVert\mathbfcal{B}\rVert_{1}+\tfrac{\tau_{B}}{2}\lVert\mathbfcal{B}-\hat{\mathbfcal{B}}^{(n)}\rVert_{\textnormal{F}}^{2}$
s.to $\displaystyle{}\
\mathbf{1}_{N_{\mathit{l}}}^{{\mathsf{H}}}\mathbf{B}_{m}=\mathbf{1}_{I_{N}}^{{\mathsf{H}}}\,,\forall
m\in\\{1,\ldots,M\\}\,,$ (14c)

where the user-defined $\tau_{X},\tau_{D},\tau_{B}\in\mathbb{R}_{++}$. Sub-
task (14c) is a composite convex minimization task under affine constraints,
and can be thus solved iteratively by [79], while (14a) and (14b) have closed-
form solutions (see Appendix C (supplementary file)).

###  IV-C Experimental setting

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Performance of MultiL-KRIM vs. other methods on different datasets
(D1, D2, D3) and sampling patterns (P1, P2). TABLE I: Average error metrics
over sampling ratios in all datasets for sampling schemes P1 and P2. The best
and second-best performing methods on each category are shown in red and blue,
respectively.

Method | Random Sampling (P1)
---|---
Sea Surface Temperature | PM 2.5 Concentration | Sea Level Pressure
MAE | RMSE | MAPE | MAE | RMSE | MAPE | MAE | RMSE | MAPE
GraphTRSS | 0.2699 | 0.4404 | 0.0250 | 2.0212 | 3.5152 | 0.3449 | 291.7323 | 494.7112 | 0.0029
KRG | 0.2097 | 0.3301 | 0.0195 | 1.7488 | 3.4578 | 0.3042 | 280.1288 | 472.0558 | 0.0027
KGL | 0.2075 | 0.3287 | 0.0188 | 1.7430 | 3.4657 | 0.3025 | 277.5922 | 461.9560 | 0.0025
LRDS | 0.1586 | 0.2559 | 0.0162 | 1.8688 | 3.4736 | 0.3438 | 271.6200 | 444.6673 | 0.0027
MMF | 0.1549 | 0.2417 | 0.0147 | 1.6751 | 3.2526 | 0.2809 | 243.4078 | 413.8005 | 0.0024
NBP | 0.1469 | 0.2242 | 0.0145 | 1.6222 | 3.1857 | 0.2726 | 241.5242 | 412.2467 | 0.0024
MultiL-KRIM | 0.1460 | 0.2315 | 0.0148 | 1.5812 | 3.0745 | 0.2743 | 231.4998 | 395.8148 | 0.0023
Method | Entire Snapshots Sampling (P2)
Sea Surface Temperature | PM 2.5 Concentration | Sea Level Pressure
MAE | RMSE | MAPE | MAE | RMSE | MAPE | MAE | RMSE | MAPE
GraphTRSS | 0.7204 | 1.1669 | 0.0751 | 2.3844 | 4.4152 | 0.4072 | 340.3902 | 550.2937 | 0.0034
KRG | 0.6744 | 1.0759 | 0.0695 | 2.1357 | 4.0605 | 0.3862 | 320.8218 | 563.9427 | 0.0032
KGL | 0.6682 | 1.0727 | 0.0691 | 2.1284 | 4.0418 | 0.3631 | 319.7752 | 553.7329 | 0.0032
LRDS | 0.6739 | 1.0754 | 0.0702 | 2.2765 | 4.1946 | 0.3724 | 331.0909 | 540.7169 | 0.0034
MMF | 0.6299 | 0.9938 | 0.0682 | 2.1077 | 4.0294 | 0.3668 | 314.2487 | 540.5447 | 0.0032
NBP | 0.6214 | 0.9875 | 0.0665 | 2.0875 | 3.9952 | 0.3756 | 310.8744 | 528.2204 | 0.0031
MultiL-KRIM | 0.5934 | 0.9566 | 0.0651 | 2.0456 | 3.8761 | 0.3738 | 306.3696 | 528.7843 | 0.0031

MultiL-KRIM is compared against several state-of-the-art methods on TVGS
recovery, such as the low-rank-and-differential-smoothness (LRDS) method [70],
the Sobolev-smoothness (GraphTRSS) one [77], the nonparametric kernel-based
KRG [55], KGL [54], and NBP [30], as well as the multi-layer matrix
factorization (MMF) [14]. Among these, LRDS, MMF, NBP, and the proposed
MultiL-KRIM promote low-rank models, while the others do not. In another view,
MultiL-KRIM, NBP, MMF, KGL, and KRG are factorization/structured methods,
while LRDS and GraphTRSS are unstructured. More readable categorization of
methods can be found at Table V (supplementary file). Notice that numerical
tests of other methods on the same datasets can be found in [69, 70, 55, 77].
MultiL-KRIM, NBP, KGL, KRG, and MMF are all solved by the successive-convex-
approximation framework of [78], and realized in software by the Julia
programming language [80]. The software code of LRDS and GraphTRSS can be
found available online. All tests were run on an 8-core Intel(R) i7-11700
2.50GHz CPU with 32GB RAM. The following sub-sections explore all steps of the
MultiL-KRIM pipeline, portraited in Figure 1.

Performance metrics. Root mean square error (RMSE), mean absolute error (MAE),
and mean absolute percentage error (MAPE) are defined in Section E-A
(supplementary file). Detailed error values of all methods across all datasets
and sampling patterns are displayed in Tables XII and XIII (supplementary
file). For all error metrics RMSE, MAE, and MAPE, lower values are better. All
methods are finely tuned to achieve the lowest MAE. The reported metric values
are the uniform averages of 20 independent runs.

Data (Step 1 of Figure 1). Three real-world datasets are considered [69, 70,
77].

  1. 1.

D1 (sea surface temperature): Monthly data from 1870 to 2014 are collected
with a spatial resolution of $1^{\circ}~{}\textnormal{latitude}\times
1^{\circ}~{}\textnormal{longitude}$ global grid. 500-month data of 100 random
locations on the Pacific ocean are selected. The dimensions of data
$\mathbf{Y}$ are $(I_{0},I_{N})=(100,500)$.

  2. 2.

D2 (particulate matter 2.5): The California daily mean PM2.5 concentration
data are collected from 93 observation sites over 200 days starting from
January 1, 2015. The dimensions of data $\mathbf{Y}$ are
$(I_{0},I_{N})=(93,200)$.

  3. 3.

D3 (sea level pressure): Global sea-level pressure data were gathered between
1948 and 2010, with a spatial resolution of
$2.5^{\circ}~{}\textnormal{latitude}\times
2.5^{\circ}~{}\textnormal{longitude}$, and temporal resolution of five days.
Here, data of 500 random sites over a span of 400 time units are considered.
The dimensions of data $\mathbf{Y}$ are $(I_{0},I_{N})=(500,400)$.

The graph topology is constructed by the $k$-nearest neighbors method as in
[70]. In particular, the Euclidean distance $d_{ij}$ between any pair of
sensors $i,j\in\mathcal{V}$ is calculated based on their given geographical
locations. Each node $i\in\mathcal{V}$ is connected to $k$ nearest neighbors
$\mathcal{N}_{i}\coloneqq\\{i_{1},i_{2},\ldots,i_{k}\\}\subset\mathcal{V}$ in
terms of the Euclidean distance Thus, $(i,i_{j})\in\mathcal{E},\forall
i_{j}\in\mathcal{N}_{i}$. Note that $i\notin\mathcal{N}_{i}$. Then, the
$(i,j)$th entry of the adjacency matrix $\mathbf{W}$ is given as
$w_{ij}=1/d^{2}_{ij}$ if $(i,j)\in\mathcal{E}$, and $w_{ij}=0$ otherwise.
Eventually, the graph Laplacian matrix $\mathbf{L}$ is calculated as described
in Section IV-A.

Sampling patterns (Step 1 of Figure 1). With sampling ratio
$r\in\\{0.1,0.2,0.3,0.4,0.5\\}$, the following sampling patterns are
considered.

  1. 1.

P1 (random sampling) [69]: signals of $\lceil I_{0}\cdot r\rceil$ randomly
sampled nodes at each time instant $t$ are observed, where $\lceil\cdot\rceil$
denotes the ceiling function.

  2. 2.

P2 (entire-snapshots sampling) [77]: all node signals are sampled at $\lceil
I_{N}\cdot r\rceil$ randomly chosen time instants. In other words, $\lceil
I_{N}\cdot r\rceil$ snapshots are randomly sampled.

Navigator data formation (Step 2 of Figure 1). The navigator data
$\check{\mathbf{Y}}_{\textnormal{nav}}$ of size $\nu\times N_{\text{nav}}$,
for some $\nu,N_{\text{nav}}\in\mathbb{N}_{*}$, are defined by the sampled-
and-zero-padded data $\mathscr{S}_{\Omega}(\mathbf{Y})$ via the following
configurations.

  1. 1.

Nav1 (snapshots): $\check{\mathbf{Y}}_{\textnormal{nav}}$ are defined by the
snapshots of $\mathscr{S}_{\Omega}(\mathbf{Y})$, represented by the columns of
$\mathscr{S}_{\Omega}(\mathbf{Y})$. For pattern P1,
$\check{\mathbf{Y}}_{\textnormal{nav}}\coloneqq\mathscr{S}_{\Omega}(\mathbf{Y})$,
with $\nu=I_{0}$ and $N_{\text{nav}}=I_{N}$; for pattern P2, missing snapshots
are discarded, i.e., $\nu=I_{0}$ and $N_{\text{nav}}=I_{N}\cdot r$.

  2. 2.

Nav2 (time profiles of nodes): $\check{\mathbf{Y}}_{\textnormal{nav}}$ are
defined by the node time profiles of $\mathscr{S}_{\Omega}(\mathbf{Y})$,
represented by the rows of $\mathscr{S}_{\Omega}(\mathbf{Y})$, i.e.,
$\check{\mathbf{Y}}_{\textnormal{nav}}\coloneqq\mathscr{S}_{\Omega}(\mathbf{Y})^{\intercal}$,
with $\nu=I_{N}$ and $N_{\text{nav}}=I_{0}$.

  3. 3.

Nav3 (patches): for the user-defined $\delta t\in\mathbb{N}$, $\delta
t<I_{N}/2$, and for each $i\in\mathcal{V}$, $t\in\mathbb{N}_{*}$, $\delta
t+1\leq t\leq I_{n}-\delta t$, a “patch”
$\mathbf{P}_{it}=\mathscr{S}_{\Omega}(\mathbf{Y})|_{\mathcal{N}_{i}\times\mathcal{N}_{t}}$
of size $k\times(2\delta t+1)$ is constructed as the submatrix of
$\mathscr{S}_{\Omega}(\mathbf{Y})$ corresponding to the rows and columns
indicated by the index sets $\mathcal{N}_{i}$ and $\mathcal{N}_{t}$,
respectively, where $\mathcal{N}_{i}$ is the $k$-nearest neighbors of node $i$
retrieved upon building $\mathbf{W}$, and
$\mathcal{N}_{t}\coloneqq\\{\tau\in\mathbb{N}_{*}{\mathbin{}\mid\mathbin{}}t-\delta
t\leq\tau\leq t+\delta t\\}$. Then, each column of the navigator data matrix
$\check{\mathbf{Y}}_{\textnormal{nav}}$ is a vectorized $\mathbf{P}_{it}$. The
size of $\check{\mathbf{Y}}_{\textnormal{nav}}$ is $k(2\delta t+1)\times
I_{0}(I_{N}-2\delta t)$. As post-processing, columns with all zeros of
$\check{\mathbf{Y}}_{\textnormal{nav}}$ are removed. This formation offers the
densest navigator point-cloud, but raises computational complexity for
landmark-point selection step.

  4. 4.

Nav4 (windows): by the same notations used in the description of Nav3, a
“window” centered at time instant $t$ is the submatrix
$\mathscr{S}_{\Omega}(\mathbf{Y})|_{\mathcal{V}\times\mathcal{N}_{t}}$. Each
column of the navigator data matrix is one vectorized window, so the size of
$\check{\mathbf{Y}}_{\textnormal{nav}}$ is $\nu=I_{0}\cdot(2\delta t+1)$ and
$N_{\text{nav}}=I_{N}-2\delta t$. Finally, drops all-zeros columns of
$\check{\mathbf{Y}}_{\textnormal{nav}}$.

Landmark point selection (Step 3 of Figure 1). Based on Euclidean distance
between navigator data points, the following landmark-point selection
strategies are examined:

  1. 1.

L1: greedy max-min-distance (maxmin) [81],

  2. 2.

L2: $k$-means clustering, and

  3. 3.

L3: fuzzy $c$-means clustering, whose centroids are chosen to be the landmark
points [82].

Kernel matrix construction (Step 4 of Figure 1). In the single-kernel case
($M=1$), the kernel from Table IV (supplementary file) which yields optimal
performance is used. The intercept $c$ of the polynomial kernel (Appendix A)
is set to be the entry-wise mean of the landmark points. After extensive
experiments, multiple kernels ($M>1$) did not score considerable improvement
over the previous single-kernel strategy.

The next sub-sections dive into different scenarios, each of which contains
four components, closely related to the steps of the pipeline. A scenario is
denoted by a tuple ($\cdot,\cdot,\cdot,\cdot$), corresponding to data,
sampling pattern, navigator data formation, and landmark-point selection
strategy. If an entry of the tuple is left empty, then all possible values of
the component are validated. For example, (D1, P1, Nav1,$\cdot$) considers sea
surface temperature, random sampling, snapshots, and validates all strategies
L1, L2, and L3.

###  IV-D Numerical tests

####  IV-D1 The sea-surface-temperature dataset

To examine the effect of different navigator data formations and landmark-
point selection strategies, tests are conducted on the (D1,
P1,$\cdot$,$\cdot$) scenario, which are detailed in Section E-B. It is
observed that Nav1 has the lowest average errors among all navigator data
formations. On the other hand, there is no significant performance difference
between landmark-point selection strategies L1, L2, and L3. Henceforth,
consider scenario ($\cdot$, P1, Nav1, L1); consider scenarios ($\cdot$, P2,
Nav1, L1) if $\lceil I_{N}\cdot r\rceil\geq 100$, and ($\cdot$, P2, Nav3, L1)
if $\lceil I_{N}\cdot r\rceil<100$.

Table I, LABEL:fig:plot.d1p1 and LABEL:fig:plot.d1p2 show that low-rank
methods outperform others under both sampling patterns. Averaged over sampling
ratios, MultiL-KRIM has lower MAE values than other methods by $1$% to $46$%
for pattern P1 and $5$% to $18$% for pattern P2. Among low-rank methods, the
kernel-based NBP and MultiL-KRIM have slight advantages over LRDS and MMF,
especially at the lowest sampling ratio of $10$% (LABEL:fig:plot.d1p1 and
LABEL:fig:plot.d1p2). The better performance of KGL and KRG over GraphTRSS
highlights the usefulness of highly structured data modeling. For pattern P2,
the gaps between all methods become narrow; this can also be seen from
LABEL:fig:plot.d1p2. Moreover, at higher sampling ratios of $40$% and $50$%,
the performance of factorization methods converges.

####  IV-D2 The PM2.5-concentration dataset

LABEL:fig:plot.d2p1 and LABEL:fig:plot.d2p2 show that unstructured methods
LRDS and GraphTRSS score the highest errors across all sampling rates in both
patterns. Unlike in the previous dataset, all structured methods outweighs
LRDS under all sampling patterns and ratios. Furthermore, the merits of
structured methods are emphasized at the lowest sampling ratio of 10%, where
they outperform LRDS and GraphTRSS by noticeable margins. Also taken from
LABEL:fig:plot.d2p1 and LABEL:fig:plot.d2p2, MultiL-KRIM, NBP, and MMF
consistently outperform KGL and KRG in both patterns, suggesting the necessity
of low-rank factorization. Compared to pattern P1, variance between all
methods is smaller in P2, especially among factorization methods, which can
also be observed from Table I. In addition, the average MAE value of MultiL-
KRIM is lower by $3$% to $22$% in pattern P1, but the range of improvement is
$3$% to $12$% in pattern P2.

####  IV-D3 The sea-level-pressure dataset

Figure 6 (supplementary file) depicts the sensitivity of MAE values to
variations of the inner dimension $d_{1}\in\\{10,20,30,40,50,60,70\\}$ and the
choice of kernel from the dictionary in Table IV (supplementary file), under
scenario (D3, P1, Nav1, L1) at sampling ratio of $10$%. The lowest MAE is
reached at $d_{1}=10$ and for the Gaussian kernel, almost regardless of its
parameter $\sigma$. Notice that the linear kernel gives the highest MAE across
all different $d_{1}$.

Table I, along with LABEL:fig:plot.d3p1 and LABEL:fig:plot.d3p2, exhibits the
continued dominance of low-rank factorization methods MultiL- KRIM, NBP, and
MMF. In terms of average MAE values, MultiL-KRIM outperforms the all competing
methods by $4$% to $21$% for pattern P1, and by $1$% to $10$% for pattern P2.
In the case of sampling pattern P1, LABEL:fig:plot.d3p1 illustrates a distinct
gap between the low-rank factorization methods MultiL-KRIM, NBP, and MMF,
compared to others. This suggests the need for both low-rank and structured
modeling. From Table I, the average MAE of LRDS in the P1 case is 3% and 4%
lower than that of KGL and KRG, respectively. This indicates a slight
advantage of low-rank modeling over kernel modeling in pattern P1. On the flip
side, LABEL:fig:plot.d3p2 shows the outperformance of structured methods
MultiL-KRIM, NBP, MMF, KGL, and KRG, over the unstructured ones LRDS and
GraphTRSS under pattern P2, advocating the interpolation strength of
structured modeling.

Computationally, MultiL-KRIM is more efficient than LRDS, since LRDS requires
SVD computations for nuclear-norm minimization, which are expected to inflict
computational burdens if the dimensions of the dataset are large. Indeed, for
dataset D3 of dimensions $500$-by-$400$, LRDS takes on average $1384$ seconds
for pattern P1 and $1456$ seconds for pattern P2, around three times the run
time of MultiL-KRIM ($444$s for P1 and $492$s for P2). Details on
computational times can be found in Table VIII (supplementary file).

##  V Application: Dynamic Magnetic Resonance Imaging

(a)

(b)

(c)

(d)

Figure 4: (a) Domain of size $I_{1}\times I_{2}\times I_{3}$ where the dMRI
data are observed and collected, with $I_{3}$ denoting the number of time
frames. The marked $\upsilon\times I_{2}\times I_{3}$ box shows the typical
location of the faithful “navigator/pilot” data, corresponding, usually, to
the “low-frequency” area of the domain. The “white dots” indicate the
locations where data are collected, with the black-colored majority of the
(k,t)-space to correspond to locations where there are missing data. After
“flattening” the 3D (k,t)-space into a 2D one (see Section II), these “white
dots” define the index set $\Omega$, and thus the sampling operator
$\mathscr{S}_{\Omega}$. (b) The $(I_{1}\times I_{2}\times I_{3})$-sized image
domain. The image-domain data are obtained by applying the two dimensional
(2D) inverse DFT $\mathscr{F}^{-1}$ to the (k,t)-space data [5]. (c) 1D
Cartesian and (d) radial sampling trajectories in k-space. The white-colored
lines indicate locations where data are collected.

This section focuses on dMRI, serving as a more detailed exposition of [53].
Section V-A defines the problem of highly-accelerated dMRI recovery, while
MultiL-KRIM’s inverse problem is stated in Section V-B. The numerical tests in
Section V-C demonstrate that MultiL-KRIM offers significant speed-ups from its
predecessor KRIM [48], with no performance losses, while outperforming state-
of-the-art “shallow” learning methods.

###  V-A DMRI-data recovery

DMRI is a common non-invasive imaging technique for tracking the motion of
body organs, offering valuable insights in the fields of cardiac and
neurological diagnosis [5]. All typical data-analytic bottlenecks show up in
dMRI: high dimensionality of the image data, many missing data entries due to
severe under-sampling to achieve high acceleration in data collection, and
strong latent spatio-temporal correlations within the data because of body-
organ structured movements; e.g., a beating heart [4, 5].

Naturally, there is a rich literature of imputation methods to address the
highly under-sampled dMRI problem: compressed sensing [83, 84], low-rank
modeling [85, 86, 87], and learning methods based on dictionaries [88, 89,
90], manifolds [91, 92, 42, 93, 94, 48], kernels [95, 94, 96, 97, 48], tensors
[98, 99], DeepL [100], and DIP [27, 101]. Several of those methods are
selected for comparisons against MultiL-KRIM in Section V-C, including the
manifold-learning-based SToRM [42], the low-rank factorization with sparsity
PS-Sparse [86], tensor factorization with total variation LRTC-TV [99], and
the DIP-based model TDDIP [27]. SToRM follows a typical ManL route and
utilizes as a regularizer a quadratic loss involving the graph Laplacian
matrix. PS-Sparse uses a “semi-blind” factorization of two orthogonal low-rank
matrices, where one of those matrix factors is built from navigator data.
LRTC-TV instead enforces a “blind” decomposition of the data matrix, and
imposes low-rankness by nuclear-norm minimization. Meanwhile, TDDIP considers
randomly generated points on a user-defined manifold (a helix) as the
navigator data.

The “(k,t)-space/domain” dMRI data are represented by a 3-way tensor
$\mathbfscr{Y}\in\mathbb{C}^{I_{1}\times I_{2}\times I_{3}}$
(LABEL:fig:dmri.ktspace), where $I_{1}$ and $I_{2}$ are the number of phase
and frequency encodings of the k-space, respectively, while $I_{3}$ stands for
the number of time frames. Following the general data formation of Section II,
$\mathbf{Y}$ is obtained by vectorizing each “slice/frame” $\mathbfcal{Y}_{t}$
of $\mathbfscr{Y}$ ($t$ denotes discrete time with
$t\in\\{1,\ldots,I_{3}\\}$), or,
$\mathbf{Y}\coloneqq[\mathbf{y}_{1},\ldots,\mathbf{y}_{I_{3}}]\in\mathbb{C}^{I_{0}\times
I_{3}}$, where $I_{0}\coloneqq I_{1}I_{2}$ and
$\mathbf{y}_{t}\coloneqq\text{vec}(\mathbfcal{Y}_{t})$. Usually, the
(k,t)-space data are severely under-sampled to promote highly accelerated data
collection [4], and to this end, popular sampling strategies are 1-D Cartesian
(LABEL:fig:cartesian.sampling) and radial (LABEL:fig:radial.sampling)
sampling, which define the index set $\Omega$ and the sampling operator
$\mathscr{S}_{\Omega}$ (Section II). The goal is to recover high fidelity
images in the image domain, despite the large number of missing entries in the
(k,t)-space. Since the end product of the dMRI recovery is an image, the
linear transform mapping $\mathscr{T}$, introduced in Section II to map data
from the image to the (k,t)-domain, becomes the 2D DFT $\mathscr{F}$.

###  V-B The inverse problem

With $\mathbf{X}$ denoting the image-domain data, the generic MultiL-KRIM’s
inverse problem (11) takes here the following special form:

$\displaystyle\min_{(\mathbf{X},\mathbf{Z},\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q},\mathbfcal{B})}{}\
$
$\displaystyle{}\tfrac{1}{2}\lVert\mathbf{X}-\mathbfcal{D}_{1}\mathbfcal{D}_{2}\cdots\mathbfcal{D}_{Q}\mathbfcal{K}\mathbfcal{B}\rVert^{2}_{\textnormal{F}}+\underbrace{\lambda_{1}\lVert\mathbfcal{B}\rVert_{1}}_{\mathcal{R}_{1}(\mathbfcal{B})}$
$\displaystyle+\underbrace{\tfrac{\lambda_{2}}{2}\lVert\mathbf{Z}-\mathscr{F}_{\text{t}}(\mathbf{X})\rVert_{\textnormal{F}}^{2}+\lambda_{3}\lVert\mathbf{Z}\rVert_{1}}_{\mathcal{R}_{21}(\mathbf{X},\mathbf{Z})}$
$\displaystyle+\underbrace{\tfrac{\lambda_{4}}{2}\sum\nolimits_{q=1}^{Q}\lVert\mathbfcal{D}_{q}\rVert_{\textnormal{F}}^{2}}_{\mathcal{R}_{22}(\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q})}$
(15a) s.to
$\displaystyle\mathscr{S}_{\Omega}(\mathbf{Y})=\mathscr{S}_{\Omega}\mathscr{F}(\mathbf{X})\,,$
(15b)
$\displaystyle\mathbf{1}_{N_{\mathit{l}}}^{{\mathsf{H}}}\mathbf{B}_{m}=\mathbf{1}_{I_{N}}^{{\mathsf{H}}}\,,\forall
m\in\\{1,\ldots,M\\}\,,$ (15c)

where the generic regularizer in (11) becomes

$\displaystyle\mathcal{R}(\mathbf{X},\mathbf{Z},\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q},\mathbfcal{B}){}\coloneqq{}$
$\displaystyle\mathcal{R}_{1}(\mathbfcal{B})+\mathcal{R}_{21}(\mathbf{X},\mathbf{Z})$
$\displaystyle+\mathcal{R}_{22}(\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q})\,.$

Regularizer $\mathcal{R}_{21}(\cdot,\cdot)$ is introduced to capitalize on the
prior knowledge that the dataset at hand, to be used in Section V-C, records
the periodic movement of a beating heart over a static background. The mapping
$\mathscr{F}_{\textnormal{t}}(\mathbf{X})$ stands for the temporal 1D DFT,
which operates on rows of $\mathbf{X}$. Each row of $\mathbf{X}$, of length
$I_{3}$, is the time profile of a single pixel in the image domain. The
periodic or static nature of each row of $\mathbf{X}$ suggests that its 1D DFT
is a sparse vector over the frequency domain. This observation justifies the
penalty term $\lVert\mathbf{Z}\rVert_{1}$, well known for promoting sparsity,
with $\mathbf{Z}$ being an auxiliary matrix-valued variable used to simplify
computations when solving the inverse problem. The solutions of the convex
sub-tasks of the inverse problem follow a similar routine to that of Section
IV-B, and are thus deferred to Appendix D (supplementary file) due to space
limitations.

###  V-C Numerical tests

(a)

(b)

Figure 5: Normalized RMSE values (the lower the better) w.r.t. the ground
truth vs. (sampling-time) acceleration factors. An acceleration factor
“$\beta$x” means that only $1/\beta$ of the full-scan $I_{1}\times I_{2}$
k-space is sampled (see LABEL:fig:cartesian.sampling and
LABEL:fig:radial.sampling).

Following [86, 93, 94, 47, 48], the proposed framework is validated on the
magnetic resonance extended cardiac-torso (MRXCAT) cine phantom dataset [102],
of size $(I_{1},I_{2},I_{3})=(408,408,360)$, under both radial and Cartesian
sampling. Similarly to [47, 48], the navigator data are formed by vectorizing
the heavily sampled region of the (k,t)-space. In other words, the
$\upsilon\times I_{2}$ box of each frame in LABEL:fig:dmri.ktspace, with
$\upsilon=4$, generates a number $I_{3}$ of $(\upsilon I_{2})\times 1$
navigator vectors. Then, the maxmin strategy [81] is adopted to select
$N_{\mathit{l}}\coloneqq 70$ landmark points. Parameter $Q\in\\{2,3\\}$, and
for each $Q$ the inner matrix dimensions of
$\\{\mathbfcal{D}_{q}\\}_{q=1}^{Q}$ are set as follows:  (i) if $Q=2$, then
$d_{1}=8$; and  (ii) if $Q=3$, then $(d_{1},d_{2})=(2,8)$.  Table IX
(supplementary file) shows the number of parameter unknowns (9) for different
settings of $Q$ to emphasize the efficiency of MultiL-KRIM’s factorization.
The inner-dimension parameter of the predecessor KRIM [48] is $d\coloneqq 8$.
Parameter $M\in\\{1,7\\}$ for both MultiL-KRIM and KRIM, with choices of
kernels as in [48] (see also Table IV (supplementary file)).

Besides the predecessors KRIM [48] and BiLMDM [47], MultiL-KRIM[$M,Q$]
competes also against the low-rank tensor factorization with total variation
(LRTC-TV) [99], the DIP-based TDDIP [27] (designed especially for radial
sampling), the popular PS-Sparse [86] and ManL-based SToRM [42]. Tags KRIM[S]
and KRIM[M] refer to [48] for the case of a single ($M=1$) and multiple
($M>1$) kernels, respectively. Performance of several other state-of-the-art
methods and their comparison against KRIM and BiLMDM on the same dataset can
be found in [47, 48]. Additionally, [48, 47] include also the validation of
KRIM and BiLMDM on several other dMRI datasets. The software code for
(MultiL-)KRIM and BiLMDM was written in Julia [80], while the software codes
of the other methods were found available online. All tests were run on an
8-core Intel(R) i7-11700 2.50GHz CPU with 32GB RAM. As mentioned after (2),
this CPU was not able to carry through NBP computations [30], since the size
$I_{0}\times I_{0}$ ($>10^{10}$ entries) of $\mathbf{K}_{\mathcal{Z}}$ seemed
to be too large for the CPU at hand.

Performance metrics. The main evaluation metric is the normalized root mean
square error
$\text{NRMSE}\coloneqq\lVert\mathbf{X}-\hat{\mathbf{X}}\rVert_{\textnormal{F}}/\lVert\mathbf{X}\rVert_{\textnormal{F}}$,
where $\mathbf{X}$ denotes the image data obtained from fully sampled
(k,t)-space data, and $\hat{\mathbf{X}}$ represents the estimate of
$\mathbf{X}$. Additionally, reconstruction quality of high-frequency regions
are quantified by the high-frequency error norm (HFEN) [88], and two sharpness
measures M1 (intensity-variance based) and M2 (energy of the image gradient)
[103, (43) and (46)]. Lastly, the structural similarity measure (SSIM) [104]
captures local similarities between the reconstructed and the ground-truth
images. Desirable outcomes are low values in NRMSE and HFEN, and high values
in SSIM, M1, and M2. All methods were finely tuned to achieve the lowest
NRMSE. All reported metric values are mean values of 10 independent runs.

It can be seen in Figure 5, and more precisely in Tables II and III, that
MultiL-KRIM matches the performance of KRIM. Nevertheless, MultiL-KRIM shows
remarkable reduction in computational time by up to $65$% for multiple kernels
(red-colored) and $30$% for a single kernel (blue-colored). To showcase that
the time reduction is because of the novel factorization approach, the
reported times of KRIM[M] and KRIM[S] do not include the time of KRIM’s
dimensionality-reduction pre-step. Furthermore, looking within MultiL-
KRIM[$M,Q$], experiments verify that compared to when $Q=2$, there is notable
speed-up when $Q=3$ and the inner dimensions $\\{d_{q}\\}_{q=1}^{Q-1}$ are
chosen appropriately. Due to intensive fine-tuning, there is only slight
difference between the single-kernel case ($M=1$) versus the multiple-kernel
case ($M=7$). The tests also demonstrate the better performance of MultiL-KRIM
over LRTC-TV, which deteriorates quickly with high acceleration rates under
Cartesian sampling, and has high computational footprint due to nuclear-norm
minimization. PS-Sparse is a contender with a trade-off between NRMSE and
runtime, but it declines noticeably at high acceleration rates, while MultiL-
KRIM does not. Additionally, structural quality measures such as SSIM and HFEN
accentuate the difference between MultiL-KRIM and PS-Sparse. Although MultiL-
KRIM does not reach the NRMSE and SSIM of TDDIP, it produces lower HFEN and
higher sharpness measure M1, i.e., sharper images. More importantly, MultiL-
KRIM offers a less computationally expensive and more data-driven approach
than TDDIP, which requires a carefully handcrafted prior (manifold structure,
here a helix) that does not explicitly utilize the observed data. Extra
experiments on TDDIP, described in Appendix F (supplementary file), show that
when the complexity of TDDIP decreases, its NRMSE and SSIM approach those of
MultiL-KRIM. Albeit, increasing the number of unknowns of MultiL-KRIM, either
with more kernels $M$ or larger inner dimensions $\\{d_{q}\\}_{q=1}^{Q-1}$,
does not improve its performance. These observations intrigue with a “deep”
architecture for MultiL-KRIM, such as adding non-linear activation function
layers into the factorization. This idea is based on the results that deep
matrix factorization outperforms the standard MMF in matrix completion tasks
[15, 16].

TABLE II: Performance on Cartesian sampling (acceleration rate: 20x)

Methods $\setminus$ Metrics | NRMSE | SSIM | HFEN | M1 | M2 | Time
---|---|---|---|---|---|---
tableMulti MultiL-KRIM[$M=7,Q=2$] | ${\mathbf{0.0443}}$ | $\mathbf{0.8698}$ | ${0.1147}$ | $37.26$ | ${1.4\times 10^{6}}$ | 1.2hrs
tableSingle MultiL-KRIM[$M=1,Q=2$] | ${0.0453}$ | ${0.8686}$ | ${0.1269}$ | $37.24$ | ${1.3\times 10^{6}}$ | 52min
tableMulti MultiL-KRIM[$M=7,Q=3$] | ${{0.0444}}$ | ${0.8695}$ | ${0.1149}$ | $37.26$ | ${1.4\times 10^{6}}$ | 1.1hrs
tableSingle MultiL-KRIM[$M=1,Q=3$] | ${0.0454}$ | ${0.8684}$ | ${0.1270}$ | $37.24$ | ${1.3\times 10^{6}}$ | 48min
tableMulti KRIM[M] [48] | $\mathbf{0.0443}$ | ${0.8696}$ | $\mathbf{0.1136}$ | $\mathbf{37.34}$ | ${1.4\times 10^{6}}$ | 3.3hrs
tableSingle KRIM[S] [48] | $0.0450$ | $0.8670$ | $0.1149$ | $37.25$ | $1.3\times 10^{6}$ | 55min
BiLMDM [47] | $0.0488$ | $0.8589$ | $0.1423$ | $37.08$ | $1.3\times 10^{6}$ | 40min
LRTC-TV [99] | $0.1645$ | $0.6943$ | $0.4330$ | $31.52$ | $\mathbf{1.5\times 10^{6}}$ | 8.3hrs
SToRM [42] | 0.0850 | $0.8110$ | $0.2504$ | $37.02$ | $1.2\times 10^{6}$ | 58min
PS-Sparse [86] | $0.0550$ | $0.8198$ | $0.1483$ | $37.18$ | $1.2\times 10^{6}$ | 15min

TABLE III: Performance on radial sampling (acceleration rate: 16x)

Methods $\setminus$ Metrics | NRMSE | SSIM | HFEN | M1 | M2 | Time
---|---|---|---|---|---|---
tableMulti MultiL-KRIM[$M=7,Q=2$] | ${{0.0448}}$ | ${0.8680}$ | ${\mathbf{0.1023}}$ | ${37.35}$ | ${1.4\times 10^{6}}$ | 1.2hrs
tableSingle MultiL-KRIM[$M=1,Q=2$] | ${0.0465}$ | ${0.8618}$ | ${0.1305}$ | $37.26$ | ${1.4\times 10^{6}}$ | 44min
tableMulti MultiL-KRIM[$M=7,Q=3$] | $0.0448$ | ${0.8675}$ | ${0.1030}$ | ${37.35}$ | ${1.4\times 10^{6}}$ | 1hr
tableSingle MultiL-KRIM[$M=1,Q=3$] | ${0.0465}$ | ${0.8618}$ | ${0.1305}$ | $37.26$ | ${1.4\times 10^{6}}$ | 41min
tableMulti KRIM[M] [48] | ${0.0450}$ | ${0.8670}$ | ${0.1136}$ | $37.34$ | ${1.4\times 10^{6}}$ | 3hrs
tableSingle KRIM[S] [48] | $0.0465$ | $0.8618$ | $0.1301$ | $37.26$ | $1.4\times 10^{6}$ | 58min
BiLMDM [47] | $0.0475$ | $0.8560$ | $0.1491$ | $37.30$ | $1.4\times 10^{6}$ | 1.4hrs
TDDIP [27] | $\mathbf{0.0376}$ | $\mathbf{0.8896}$ | $0.1452$ | ${36.24}$ | $1.3\times 10^{6}$ | 3.9hrs
LRTC-TV [99] | $0.0738$ | $0.8063$ | $0.3725$ | $37.06$ | $\mathbf{1.7\times 10^{6}}$ | 8.6hrs
SToRM [42] | $0.0753$ | $0.8319$ | $0.3694$ | $\mathbf{37.38}$ | $1.6\times 10^{6}$ | 30min
PS-Sparse [86] | $0.0496$ | $0.7908$ | $0.1733$ | $37.31$ | $1.4\times 10^{6}$ | 15min

##  VI Conclusions

This paper introduced MultiL-KRIM, a nonparametric (kernel-based) data-
imputation-by-regression framework. Two important application domains were
considered for validation: time-varying-graph-signal (TVGS) recovery, and
reconstruction of highly accelerated dynamic-magnetic-resonance-imaging (dMRI)
data. Numerical tests on TVGS recovery demonstrated MultiL-KRIM’s
outperformance over state-of-the-art “shallow” methods. Furthermore, via tests
on dMRI data, MultiL-KRIM was shown to match the performance of its
predecessor KRIM [48] with significant reduction in computational time.
MultiL-KRIM outperformed state-of-the-art “shallow” methods also on dMRI data,
and although it did not reach the accuracy of the deep-image-prior-based TDDIP
[27], it offered a more explainable and computationally tractable path.
Prospectively, to further enhance MultiL-KRIM’s performance, current research
effort includes also its “deep” variations.

## References

  * [1] A. Ramlatchan, M. Yang, Q. Liu, M. Li, J. Wang, and Y. Li, “A survey of matrix completion methods for recommendation systems,” Big Data Mining and Analytics, vol. 1, no. 4, pp. 308–323, 2018.
  * [2] H. Shen, X. Li, Q. Cheng, C. Zeng, G. Yang, H. Li, and L. Zhang, “Missing information reconstruction of remote sensing data: A technical review,” IEEE Geoscience and Remote Sensing Magazine, vol. 3, no. 3, pp. 61–85, 2015\.
  * [3] J. Du, M. Hu, and W. Zhang, “Missing data problem in the monitoring system: A review,” IEEE Sensors Journal, vol. 20, no. 23, pp. 13984–13998, 2020\.
  * [4] Z.-P. Liang and P. C. Lauterbur, “An efficient method for dynamic magnetic resonance imaging,” IEEE Trans. Medical Imaging, vol. 13, no. 4, pp. 677–686, 1994.
  * [5] Z.-P. Liang and P. C. Lauterbur, Principles of Magnetic Resonance Imaging: A Signal Processing Perspective.  IEEE Press, 2000.
  * [6] M. S. Osman, A. M. Abu-Mahfouz, and P. R. Page, “A survey on data imputation techniques: Water distribution system as a use case,” IEEE Access, vol. 6, pp. 63279–63291, 2018.
  * [7] E. J. Candes and B. Recht, “Exact low-rank matrix completion via convex optimization,” in 2008 46th Annual Allerton Conference on Communication, Control, and Computing, pp. 806–812, IEEE, 2008.
  * [8] P. Chen and D. Suter, “Recovering the missing components in a large noisy low-rank matrix: Application to SFM,” IEEE transactions on pattern analysis and machine intelligence, vol. 26, no. 8, pp. 1051–1063, 2004.
  * [9] Y. Koren, R. Bell, and C. Volinsky, “Matrix factorization techniques for recommender systems,” Computer, vol. 42, no. 8, pp. 30–37, 2009.
  * [10] J. Mairal, J. Ponce, G. Sapiro, A. Zisserman, and F. Bach, “Supervised dictionary learning,” Advances in neural information processing systems, vol. 21, 2008.
  * [11] I. Tošić and P. Frossard, “Dictionary learning,” IEEE Signal Processing Magazine, vol. 28, no. 2, pp. 27–38, 2011.
  * [12] J. Sun, Q. Qu, and J. Wright, “Complete dictionary recovery over the sphere I: Overview and the geometric picture,” IEEE Transactions on Information Theory, vol. 63, no. 2, pp. 853–884, 2016.
  * [13] J. Sun, Q. Qu, and J. Wright, “Complete dictionary recovery over the sphere II: Recovery by riemannian trust-region method,” IEEE Transactions on Information Theory, vol. 63, no. 2, pp. 885–914, 2016.
  * [14] A. Cichocki and R. Zdunek, “Multilayer nonnegative matrix factorization using projected gradient approaches,” International Journal of Neural Systems, vol. 17, no. 06, pp. 431–446, 2007.
  * [15] J. Fan and J. Cheng, “Matrix completion by deep matrix factorization,” Neural Networks, vol. 98, pp. 34–41, 2018.
  * [16] P. De Handschutter, N. Gillis, and X. Siebert, “A survey on deep matrix factorizations,” Computer Science Review, vol. 42, p. 100423, 2021.
  * [17] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.
  * [18] I. V. Oseledets, “Tensor-train decomposition,” SIAM Journal on Scientific Computing, vol. 33, no. 5, pp. 2295–2317, 2011.
  * [19] Q. Zhao, G. Zhou, S. Xie, L. Zhang, and A. Cichocki, “Tensor ring decomposition,” arXiv preprint arXiv:1606.05535, 2016.
  * [20] Y.-B. Zheng, T.-Z. Huang, X.-L. Zhao, Q. Zhao, and T.-X. Jiang, “Fully-connected tensor network decomposition and its application to higher-order tensor completion,” in Proceedings of the AAAI conference on artificial intelligence, vol. 35, pp. 11071–11078, 2021.
  * [21] L. T. Nguyen, J. Kim, and B. Shim, “Low-rank matrix completion: A contemporary survey,” IEEE Access, vol. 7, pp. 94215–94237, 2019.
  * [22] J. Su, B. Xu, and H. Yin, “A survey of deep learning approaches to image restoration,” Neurocomputing, vol. 487, pp. 46–65, 2022.
  * [23] V. Antun, F. Renna, C. Poon, B. Adcock, and A. C. Hansen, “On instabilities of deep learning in image reconstruction and the potential costs of AI,” Proceedings of the National Academy of Sciences, vol. 117, no. 48, pp. 30088–30095, 2020.
  * [24] D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Deep image prior,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 9446–9454, 2018.
  * [25] R. Heckel and M. Soltanolkotabi, “Compressive sensing with un-trained neural networks: Gradient descent finds a smooth approximation,” in International Conference on Machine Learning, pp. 4149–4158, PMLR, 2020.
  * [26] S. Rey, S. Segarra, R. Heckel, and A. G. Marques, “Untrained graph neural networks for denoising,” IEEE Transactions on Signal Processing, vol. 70, pp. 5708–5723, 2022.
  * [27] J. Yoo, K. H. Jin, H. Gupta, J. Yerly, M. Stuber, and M. Unser, “Time-dependent deep image prior for dynamic MRI,” IEEE Trans. Medical Imaging, vol. 40, no. 12, pp. 3337–3348, 2021.
  * [28] Q. Zou, A. H. Ahmed, P. Nagpal, S. Kruger, and M. Jacob, “Dynamic imaging using a deep generative SToRM (Gen-SToRM) model,” IEEE Trans. Medical Imaging, vol. 40, no. 11, pp. 3102–3112, 2021.
  * [29] L. Györfi, M. Kohler, A. Krzyżak, and H. Walk, A Distribution-Free Theory of Nonparametric Regression.  New York: Springer, 2010.
  * [30] J. A. Bazerque and G. B. Giannakis, “Nonparametric basis pursuit via sparse kernel-based learning: A unifying view with advances in blind methods,” IEEE Signal Processing Magazine, vol. 30, no. 4, pp. 112–125, 2013.
  * [31] J. Fan and T. W. Chow, “Non-linear matrix completion,” Pattern Recognition, vol. 77, pp. 378–394, 2018.
  * [32] M. Stock, T. Pahikkala, A. Airola, B. De Baets, and W. Waegeman, “A comparative study of pairwise learning methods based on kernel ridge regression,” Neural computation, vol. 30, no. 8, pp. 2245–2283, 2018.
  * [33] P. Giménez-Febrer, A. Pagès-Zamora, and G. B. Giannakis, “Matrix completion and extrapolation via kernel regression,” IEEE Transactions on Signal Processing, vol. 67, no. 19, pp. 5004–5017, 2019.
  * [34] Z. Chen, W. Zhao, and S. Wang, “Kernel meets recommender systems: A multi-kernel interpolation for matrix completion,” Expert Systems with Applications, vol. 168, p. 114436, 2021.
  * [35] N. Aronszajn, “Theory of reproducing kernels,” Trans. American Mathematical Society, vol. 68, no. 3, pp. 337–404, 1950.
  * [36] B. Schölkopf and A. J. Smola, Learning with Kernels.  MIT press, 2002.
  * [37] T. Lin and H. Zha, “Riemannian manifold learning,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 30, no. 5, pp. 796–809, 2008.
  * [38] S. T. Roweis and L. K. Saul, “Nonlinear dimensionality reduction by locally linear embedding,” Science, vol. 290, no. 5500, pp. 2323–2326, 2000.
  * [39] M. Belkin and P. Niyogi, “Laplacian eigenmaps and spectral techniques for embedding and clustering,” NIPS, vol. 14, 2001.
  * [40] J. B. Tenenbaum, V. d. Silva, and J. C. Langford, “A global geometric framework for nonlinear dimensionality reduction,” Science, vol. 290, no. 5500, pp. 2319–2323, 2000.
  * [41] V. Kalofolias, X. Bresson, M. Bronstein, and P. Vandergheynst, “Matrix completion on graphs,” arXiv preprint arXiv:1408.1717, 2014.
  * [42] S. Poddar and M. Jacob, “Dynamic MRI using smoothness regularization on manifolds (SToRM),” IEEE Trans. Medical Imaging, vol. 35, no. 4, pp. 1106–1115, 2016.
  * [43] H. Ma, D. Zhou, C. Liu, M. R. Lyu, and I. King, “Recommender systems with social regularization,” in Proceedings of the fourth ACM international conference on Web search and data mining, pp. 287–296, 2011.
  * [44] N. Rao, H.-F. Yu, P. K. Ravikumar, and I. S. Dhillon, “Collaborative filtering with graph information: Consistency and scalable methods,” Advances in neural information processing systems, vol. 28, 2015.
  * [45] A. Mongia and A. Majumdar, “Matrix completion on multiple graphs: Application in collaborative filtering,” Signal Processing, vol. 165, pp. 144–148, 2019.
  * [46] Z. Shi, S. Osher, and W. Zhu, “Weighted nonlocal Laplacian on interpolation from sparse data,” Journal of Scientific Computing, vol. 73, pp. 1164–1177, 2017.
  * [47] G. N. Shetty, K. Slavakis, A. Bose, U. Nakarmi, G. Scutari, and L. Ying, “Bi-linear modeling of data manifolds for dynamic-MRI recovery,” IEEE Trans. Medical Imaging, vol. 39, no. 3, pp. 688–702, 2020.
  * [48] K. Slavakis, G. N. Shetty, L. Cannelli, G. Scutari, U. Nakarmi, and L. Ying, “Kernel regression imputation in manifolds via bi-linear modeling: The dynamic-MRI case,” IEEE Trans. Computational Imaging, vol. 8, pp. 133–147, 2022.
  * [49] M. A. Davenport and J. Romberg, “An overview of low-rank matrix recovery from incomplete observations,” IEEE Journal of Selected Topics in Signal Processing, vol. 10, no. 4, pp. 608–622, 2016.
  * [50] Y. Chi, Y. M. Lu, and Y. Chen, “Nonconvex optimization meets low-rank matrix factorization: An overview,” IEEE Transactions on Signal Processing, vol. 67, no. 20, pp. 5239–5269, 2019.
  * [51] P. Zhou, C. Lu, Z. Lin, and C. Zhang, “Tensor factorization for low-rank tensor completion,” IEEE Transactions on Image Processing, vol. 27, no. 3, pp. 1152–1163, 2017.
  * [52] C. Cai, G. Li, H. V. Poor, and Y. Chen, “Nonconvex low-rank tensor completion from noisy data,” Advances in neural information processing systems, vol. 32, 2019.
  * [53] D. T. Nguyen and K. Slavakis, “Multi-linear kernel regression and imputation via manifold learning: The dynamic MRI case,” in Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), (Seoul: Korea), April 14–19 2024.  to be presented.
  * [54] X. Pu, S. L. Chau, X. Dong, and D. Sejdinovic, “Kernel-based graph learning from smooth signals: A functional viewpoint,” IEEE Transactions on Signal and Information Processing over Networks, vol. 7, pp. 192–207, 2021.
  * [55] A. Venkitaraman, S. Chatterjee, and P. Händel, “Predicting graph signals using kernel regression where the input signal is agnostic to a graph,” IEEE Transactions on Signal and Information Processing over Networks, vol. 5, no. 4, pp. 698–710, 2019.
  * [56] M. A. Carreira-Perpin, Z. Lu, et al., “Manifold learning and missing data recovery through unsupervised regression,” in Proc. IEEE International Conference on Data Mining, pp. 1014–1019, IEEE, 2011.
  * [57] J. W. Robbin and D. A. Salamon, Introduction to Differential Geometry.  Berlin: Springer, 2022.
  * [58] A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” Proceedings of the IEEE, vol. 106, no. 5, pp. 808–828, 2018\.
  * [59] W. Hu, J. Pang, X. Liu, D. Tian, C.-W. Lin, and A. Vetro, “Graph signal processing for geometric data and beyond: Theory and applications,” IEEE Transactions on Multimedia, vol. 24, pp. 3961–3977, 2021.
  * [60] W. Huang, A. G. Marques, and A. R. Ribeiro, “Rating prediction via graph signal processing,” IEEE Transactions on Signal Processing, vol. 66, no. 19, pp. 5066–5081, 2018.
  * [61] R. Ramakrishna and A. Scaglione, “Grid-graph signal processing (grid-GSP): A graph signal processing framework for the power grid,” IEEE Transactions on Signal Processing, vol. 69, pp. 2725–2739, 2021.
  * [62] I. Jabłoński, “Graph signal processing in applications to sensor networks, smart grids, and smart cities,” IEEE Sensors Journal, vol. 17, no. 23, pp. 7659–7666, 2017.
  * [63] P.-A. Laharotte, R. Billot, E. Come, L. Oukhellou, A. Nantes, and N.-E. El Faouzi, “Spatiotemporal analysis of bluetooth data: Application to a large urban network,” IEEE Transactions on Intelligent Transportation Systems, vol. 16, no. 3, pp. 1439–1448, 2014.
  * [64] B. Płaczek, “Selective data collection in vehicular networks for traffic control applications,” Transportation Research Part C: Emerging Technologies, vol. 23, pp. 14–28, 2012.
  * [65] J. Cheng, Q. Ye, H. Jiang, D. Wang, and C. Wang, “STCDG: An efficient data gathering algorithm based on matrix completion for wireless sensor networks,” IEEE Transactions on Wireless Communications, vol. 12, no. 2, pp. 850–861, 2012.
  * [66] S. Chen, A. Sandryhaila, J. M. Moura, and J. Kovačević, “Signal recovery on graphs: Variation minimization,” IEEE Transactions on Signal Processing, vol. 63, no. 17, pp. 4609–4624, 2015.
  * [67] S. Chen, R. Varma, A. Singh, and J. Kovačević, “Signal recovery on graphs: Fundamental limits of sampling strategies,” IEEE Transactions on Signal and Information Processing over Networks, vol. 2, no. 4, pp. 539–554, 2016.
  * [68] D. Romero, V. N. Ioannidis, and G. B. Giannakis, “Kernel-based reconstruction of space-time functions on dynamic graphs,” IEEE Journal of Selected Topics in Signal Processing, vol. 11, no. 6, pp. 856–869, 2017.
  * [69] K. Qiu, X. Mao, X. Shen, X. Wang, T. Li, and Y. Gu, “Time-varying graph signal reconstruction,” IEEE Journal of Selected Topics in Signal Processing, vol. 11, no. 6, pp. 870–883, 2017.
  * [70] X. Mao, K. Qiu, T. Li, and Y. Gu, “Spatio-temporal signal recovery based on low rank and differential smoothness,” IEEE Transactions on Signal Processing, vol. 66, no. 23, pp. 6281–6296, 2018.
  * [71] M. Roughan, Y. Zhang, W. Willinger, and L. Qiu, “Spatio-temporal compressive sensing and internet traffic matrices (extended version),” IEEE/ACM Transactions on Networking, vol. 20, no. 3, pp. 662–676, 2011.
  * [72] L. Kong, M. Xia, X.-Y. Liu, M.-Y. Wu, and X. Liu, “Data loss and reconstruction in sensor networks,” in 2013 Proceedings IEEE INFOCOM, pp. 1654–1662, IEEE, 2013.
  * [73] K. Xie, X. Ning, X. Wang, D. Xie, J. Cao, G. Xie, and J. Wen, “Recover corrupted data in sensor networks: A matrix completion solution,” IEEE Transactions on Mobile Computing, vol. 16, no. 5, pp. 1434–1448, 2016.
  * [74] X. Song, Y. Guo, N. Li, and S. Yang, “A novel approach based on matrix factorization for recovering missing time series sensor data,” IEEE Sensors Journal, vol. 20, no. 22, pp. 13491–13500, 2020.
  * [75] M. Lei, A. Labbe, Y. Wu, and L. Sun, “Bayesian kernelized matrix factorization for spatiotemporal traffic data imputation and kriging,” IEEE Transactions on Intelligent Transportation Systems, vol. 23, no. 10, pp. 18962–18974, 2022.
  * [76] X. Piao, Y. Hu, Y. Sun, B. Yin, and J. Gao, “Correlated spatio-temporal data collection in wireless sensor networks based on low rank matrix approximation and optimized node sampling,” Sensors, vol. 14, no. 12, pp. 23137–23158, 2014.
  * [77] J. H. Giraldo, A. Mahmood, B. Garcia-Garcia, D. Thanou, and T. Bouwmans, “Reconstruction of time-varying graph signals via Sobolev smoothness,” IEEE Transactions on Signal and Information Processing over Networks, vol. 8, pp. 201–214, 2022.
  * [78] F. Facchinei, G. Scutari, and S. Sagratella, “Parallel selective algorithms for nonconvex big data optimization,” IEEE Trans. Signal Processing, vol. 63, no. 7, pp. 1874–1889, 2015.
  * [79] K. Slavakis and I. Yamada, “Fejér-monotone hybrid steepest descent method for affinely constrained and composite convex minimization tasks,” Optimization, vol. 67, no. 11, pp. 1963–2001, 2018.
  * [80] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM Review, vol. 59, no. 1, pp. 65–98, 2017.
  * [81] V. De Silva and J. B. Tenenbaum, “Sparse multidimensional scaling using landmark points,” tech. rep., Stanford University, 2004.
  * [82] M. Rashmi and P. Sankaran, “Optimal landmark point selection using clustering for manifold modeling and data classification,” Journal of Classification, vol. 36, no. 1, pp. 94–112, 2019.
  * [83] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182–1195, 2007.
  * [84] L. Feng, L. Axel, H. Chandarana, K. T. Block, D. K. Sodickson, and R. Otazo, “XD-GRASP: Golden-angle radial MRI with reconstruction of extra motion-state dimensions using compressed sensing,” Magnetic Resonance in Medicine, vol. 75, no. 2, pp. 775–788, 2016.
  * [85] S. G. Lingala, Y. Hu, E. V. R. DiBella, and M. Jacob, “Accelerated dynamic MRI exploiting sparsity and low-rank structure: k-t SLR,” IEEE Trans. Medical Imaging, vol. 30, no. 5, pp. 1042–1054, 2011.
  * [86] B. Zhao, J. P. Haldar, A. G. Christodoulou, and Z.-P. Liang, “Image reconstruction from highly undersampled (k,t)-space data with joint partial separability and sparsity constraints,” IEEE Trans. Medical Imaging, vol. 31, no. 9, pp. 1809–1820, 2012.
  * [87] R. Otazo, E. Candes, and D. K. Sodickson, “Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components,” Magnetic Resonance in Medicine, vol. 73, no. 3, pp. 1125–1136, 2015.
  * [88] S. Ravishankar and Y. Bresler, “MR image reconstruction from highly undersampled k-space data by dictionary learning,” IEEE Trans. Medical Imaging, vol. 30, no. 5, pp. 1028–1041, 2010.
  * [89] Y. Wang and L. Ying, “Compressed sensing dynamic cardiac cine MRI using learned spatiotemporal dictionary,” IEEE Trans. Biomedical Eng., vol. 61, no. 4, pp. 1109–1120, 2014.
  * [90] J. Caballero, A. N. Price, D. Rueckert, and J. V. Hajnal, “Dictionary learning and time sparsity for dynamic MR data reconstruction,” IEEE Trans. Medical Imaging, vol. 33, no. 4, pp. 979–994, 2014.
  * [91] M. Usman, D. Atkinson, C. Kolbitsch, T. Schaeffter, and C. Prieto, “Manifold learning based ECG-free free-breathing cardiac CINE MRI,” J. Magnetic Resonance Imag., vol. 41, no. 6, pp. 1521–1527, 2015.
  * [92] X. Chen, M. Usman, C. F. Baumgartner, D. R. Balfour, P. K. Marsden, A. J. Reader, C. Prieto, and A. P. King, “High-resolution self-gated dynamic abdominal MRI using manifold alignment,” IEEE Trans. Medical Imaging, vol. 36, no. 4, pp. 960–971, 2017.
  * [93] U. Nakarmi, K. Slavakis, J. Lyu, and L. Ying, “M-MRI: A manifold-based framework to highly accelerated dynamic magnetic resonance imaging,” in Proceedings of ISBI, pp. 19–22, 2017.
  * [94] U. Nakarmi, K. Slavakis, and L. Ying, “MLS: Joint manifold-learning and sparsity-aware framework for highly accelerated dynamic magnetic resonance imaging,” in Proceedings of ISBI, pp. 1213–1216, 2018.
  * [95] U. Nakarmi, Y. Wang, J. Lyu, D. Liang, and L. Ying, “A kernel-based low-rank (KLR) model for low-dimensional manifold recovery in highly accelerated dynamic MRI,” IEEE Trans. Medical Imaging, vol. 36, no. 11, pp. 2297–2307, 2017.
  * [96] S. Poddar, Y. Q. Mohsin, D. Ansah, B. Thattaliyath, R. Ashwath, and M. Jacob, “Manifold recovery using kernel low-rank regularization: application to dynamic imaging,” IEEE Trans. Computational Imaging, vol. 5, no. 3, pp. 478–491, 2019.
  * [97] O. Arif, H. Afzal, H. Abbas, M. F. Amjad, J. Wan, and R. Nawaz, “Accelerated dynamic MRI using kernel-based low rank constraint,” J. Medical Systems, vol. 43, no. 8, p. 271, 2019.
  * [98] J. Liu, P. Musialski, P. Wonka, and J. Ye, “Tensor completion for estimating missing values in visual data,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 35, no. 1, pp. 208–220, 2012.
  * [99] X. Li, Y. Ye, and X. Xu, “Low-rank tensor completion with total variation for visual data inpainting,” Proceedings of AAAI Conference on Artificial Intelligence, vol. 31, Feb. 2017.
  * [100] F. Knoll, K. Hammernik, C. Zhang, S. Moeller, T. Pock, D. K. Sodickson, and M. Akcakaya, “Deep-learning methods for parallel magnetic resonance imaging reconstruction: A survey of the current approaches, trends, and issues,” IEEE signal processing magazine, vol. 37, no. 1, pp. 128–140, 2020.
  * [101] A. Qayyum, I. Ilahi, F. Shamshad, F. Boussaid, M. Bennamoun, and J. Qadir, “Untrained neural network priors for inverse imaging problems: A survey,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 2022.
  * [102] L. Wissmann, C. Santelli, W. P. Segars, and S. Kozerke, “MRXCAT: Realistic numerical phantoms for cardiovascular magnetic resonance,” J. Cardiovascular Magnetic Resonance, vol. 16, no. 1, pp. 1–11, 2014.
  * [103] M. Subbarao, T.-S. Choi, and A. Nikzad, “Focusing techniques,” Optical Engineering, vol. 32, no. 11, pp. 2824–2836, 1993.
  * [104] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. Medical Imaging, vol. 13, no. 4, pp. 600–612, 2004.
  * [105] I. Steinwart, D. Hush, and C. Scovel, “An explicit description of the reproducing kernel Hilbert spaces of Gaussian RBF kernels,” IEEE Transactions on Information Theory, vol. 52, no. 10, pp. 4635–4643, 2006.
  * [106] P. Bouboulis and S. Theodoridis, “Extension of Wirtinger’s calculus to reproducing kernel Hilbert spaces and the complex kernel LMS,” IEEE Trans. Signal Processing, vol. 59, no. 3, pp. 964–978, 2010.
  * [107] K. Slavakis, P. Bouboulis, and S. Theodoridis, Chapter 17: Online Learning in Reproducing Kernel Hilbert Spaces, vol. 1 of Academic Press Library in Signal Processing: Volume 1, pp. 883–987.  Elsevier, Jan. 2014.
  * [108] H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces.  New York: Springer, 2011.

Supplementary File

##  Appendix A Reproducing Kernel Hilbert Spaces

A feature mapping $\varphi:\mathbb{C}^{\nu}\to\mathscr{H}$ maps a vector
$\mathbfit{l}\in\mathbb{C}^{\nu}$ to $\varphi(\mathbfit{l})$ in a feature
space $\mathscr{H}$. To provide structured solutions, it is assumed that
$\mathscr{H}$ is a complex-valued reproducing kernel Hilbert space (RKHS),
associated with a reproducing kernel
$\kappa(\cdot,\cdot):\mathbb{C}^{\nu}\times\mathbb{C}^{\nu}\to\mathbb{C}$,
with rich attributes in approximation theory [35, 36, 105, 106, 107]. In the
RKHS setting, the feature mapping becomes
$\varphi(\mathbfit{l})\coloneqq\kappa(\mathbfit{l},\cdot)\in\mathscr{H}$,
$\forall\mathbfit{l}\in\mathbb{C}^{\nu}$. Most well-known kernels for RKHSs
are:  (i) The linear kernel
$\kappa_{\text{L}}(\mathbfit{l},\mathbfit{l}^{\prime})\coloneqq\mathbfit{l}^{\mathsf{H}}\mathbfit{l}^{\prime}$,
where $\mathsf{H}$ denotes complex conjugate vector/matrix transposition;
(ii) the Gaussian kernel
$\kappa_{\text{G};\gamma}(\mathbfit{l},\mathbfit{l}^{\prime})\coloneqq\text{exp}[-\gamma(\mathbfit{l}-\overline{\mathbfit{l}}^{\prime})^{\intercal}(\mathbfit{l}-\overline{\mathbfit{l}}^{\prime})]$
[105, 106], where $\gamma\in\mathbb{R}_{>0}$, superscript $\intercal$ stands
for vector/matrix transposition, and the overline symbol denotes entry-wise
complex conju- gation of a scalar/vector/matrix; and  (iii) the polynomial
kernel
$\kappa_{\text{P};(c,r)}(\mathbfit{l},\mathbfit{l}^{\prime})\coloneqq(\mathbfit{l}^{{\mathsf{H}}}\mathbfit{l}^{\prime}+c)^{r}$,
where $r\in\mathbb{N}_{*}$, $c\in\mathbb{C}$.  The list of kernels used for
numerical tests is shown in Table IV.

TABLE IV: Dictionary of kernels.  Kernel index $m$ | Kernel function | Parameters
---|---|---
1 | Gaussian | $\sigma=0.2$
2 | Gaussian | $\sigma=0.4$
3 | Gaussian | $\sigma=0.8$
4 | Polynomial | $r=1$
5 | Polynomial | $r=2$
6 | Polynomial | $r=3$
7 | Polynomial | $r=4$

##  Appendix B Derivation of a compact MultiL-KRIM data modeling

$\displaystyle\mathbf{X}\approx\underbrace{\left[\begin{smallmatrix}\mathbf{D}_{1}^{(1)}&\mathbf{D}_{2}^{(1)}&\ldots&\mathbf{D}_{M}^{(1)}\end{smallmatrix}\right]}_{\mathbfcal{D}_{1}}\underbrace{\left[\begin{smallmatrix}\mathbf{D}_{1}^{(2)}&0&\cdots&0\\\
0&\mathbf{D}_{2}^{(2)}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\mathbf{D}_{M}^{(2)}\\\
\end{smallmatrix}\right]}_{\mathbfcal{D}_{2}}\ldots\underbrace{\left[\begin{smallmatrix}\mathbf{D}_{1}^{(Q)}&0&\cdots&0\\\
0&\mathbf{D}_{2}^{(Q)}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\mathbf{D}_{M}^{(Q)}\\\
\end{smallmatrix}\right]}_{\mathbfcal{D}_{Q}}\underbrace{\left[\begin{smallmatrix}\mathbf{K}_{1}&0&\cdots&0\\\
0&\mathbf{K}_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\mathbf{K}_{M}\\\
\end{smallmatrix}\right]}_{\mathbfcal{K}}\underbrace{\left[\begin{smallmatrix}\mathbf{B}_{1}\\\
\mathbf{B}_{2}\\\ \vdots\\\
\mathbf{B}_{M}\end{smallmatrix}\right]}_{\mathbfcal{B}}\,.$ (16)

By defining the “supermatrices”

$\displaystyle\mathbfcal{D}_{1}$
$\displaystyle{}\coloneqq{}[\mathbf{D}_{1}^{(1)},\mathbf{D}_{2}^{(1)},\ldots,\mathbf{D}_{M}^{(1)}]$
$\displaystyle\in\mathbb{C}^{I_{0}\times d_{1}M}\,,$
$\displaystyle\mathbfcal{D}_{q}$
$\displaystyle\coloneqq\operatorname{bdiag}{(\mathbf{D}_{1}^{(q)},\mathbf{D}_{2}^{(q)},\ldots,\mathbf{D}_{M}^{(q)})}$
$\displaystyle\in\mathbb{C}^{d_{q-1}M\times d_{q}M}\,,$
$\displaystyle\hphantom{{}\coloneqq{}}\forall q\in\\{2,\ldots,Q\\}\,,$
$\displaystyle\mathbfcal{K}$
$\displaystyle\coloneqq\operatorname{bdiag}{(\mathbf{K}_{1},\mathbf{K}_{2},\ldots,\mathbf{K}_{M})}$
$\displaystyle\in\mathbb{C}^{MN_{\mathit{l}}\times MN_{\mathit{l}}}\,,$ (17)
$\displaystyle\mathbfcal{B}$
$\displaystyle\coloneqq[\mathbf{B}_{1}^{\mathsf{H}},\mathbf{B}_{2}^{\mathsf{H}},\ldots,\mathbf{B}_{M}^{\mathsf{H}}]^{\mathsf{H}}$
$\displaystyle\in\mathbb{C}^{MN_{\mathit{l}}\times I_{N}}\,,$

data modeling (8) of MultiL-KRIM, incorporated in the generic inverse problem
(11), takes the form in Eq. 16.

##  Appendix C Solving the TVGS inverse problem

This section details the solutions to the TVGS sub-tasks (14a) and (14b). Sub-
task (14c) involves a composite convex loss function with affine constraints
which can be efficiently solved by [79].

###  C-A Solving (14a)

Given the index set $\Omega$, define the affine constraint
$\mathcal{A}_{\Omega}\coloneqq\\{\mathbf{X}\in\mathbb{C}^{I_{0}\times
I_{N}}{\mathbin{}\mid\mathbin{}}\mathscr{S}_{\Omega}(\mathbf{Y})=\mathscr{S}_{\Omega}(\mathbf{X})\\}$,
and notice that (14a) takes the form:

$\displaystyle\hat{\mathbf{X}}^{(n+1/2)}=\arg\min\nolimits_{\mathbf{X}\in\mathcal{A}_{\Omega}}\mathcal{L}_{X}(\mathbf{X})\,,$
(18)

where the differentiable strongly convex loss

$\displaystyle\mathcal{L}_{X}(\mathbf{X})$ $\displaystyle{}\coloneqq{}$
$\displaystyle\tfrac{1}{2}\lVert{\mathbf{X}}-\hat{\mathbf{M}}^{(n)}\rVert_{\textnormal{F}}^{2}+\tfrac{\lambda_{L}}{2}\operatorname{tr}(\mathbf{X}^{\intercal}\mathbf{L}_{\epsilon}^{\beta}\mathbf{X}\mathbf{\Delta}\mathbf{\Delta}^{\intercal})$
$\displaystyle+\tfrac{\tau_{X}}{2}\lVert\mathbf{X}-\hat{\mathbf{X}}^{(n)}\rVert_{\textnormal{F}}^{2}\,,$

and
$\hat{\mathbf{M}}^{(n)}\coloneqq\hat{\mathbfcal{D}}_{1}^{(n)}\cdots\hat{\mathbfcal{D}}_{Q}^{(n)}\mathbfcal{K}\hat{\mathbfcal{B}}^{(n)}$.
For convenience in notations, let
$\mathbf{X}_{*}\coloneqq\hat{\mathbf{X}}^{(n+1/2)}$.

Problem (18) is a convex minimization task with affine constraints, so it can
be solved iteratively by the efficient [79]. Nevertheless, (18) possesses also
a closed-form solution. To formulate that closed-form solution, define the
sampling matrix $\mathbf{S}_{\Omega}\in\mathbb{C}^{I_{0}\times I_{N}}$ as
follows: $[\mathbf{S}_{\Omega}]_{ij}\coloneqq 1$, if $(i,j)\in\Omega$, while
$[\mathbf{S}_{\Omega}]_{ij}\coloneqq 0$, if $(i,j)\notin\Omega$. Then,
$\mathscr{S}_{\Omega}(\mathbf{X})=\mathbf{S}_{\Omega}\odot\mathbf{X}$, where
$\odot$ stands for the Hadamard product. Furthermore, define the complement
set
$\Omega^{\complement}\coloneqq\\{(i,j)\in\\{1,\ldots,I_{0}\\}\times\\{1,\ldots,I_{N}\\}\nonscript\,|\allowbreak\nonscript\,\mathopen{}(i,j)\notin\Omega\\}$.
It can be verified that any $\mathbf{X}\in\mathbb{C}^{I_{0}\times I_{N}}$ can
be written as
$\mathbf{X}=\mathscr{S}_{\Omega}(\mathbf{X})+\mathscr{S}_{\Omega^{\complement}}(\mathbf{X})=\mathbf{S}_{\Omega}\odot\mathbf{X}+\mathbf{S}_{\Omega^{\complement}}\odot\mathbf{X}$.
Moreover, let $\operatorname{vec}\colon\mathbb{C}^{I_{0}\times
I_{N}}\to\mathbb{C}^{I_{0}I_{N}}$ denote the standard vectorization operator
which turns matrices into vectors, and let its inverse mapping be denoted by
$\operatorname{vec}^{-1}\colon\mathbb{C}^{I_{0}I_{N}}\to\mathbb{C}^{I_{0}\times
I_{N}}$. Let also the “vectorization” of $\Omega$ as the one-dimensional index
set
$\operatorname{vec}\Omega\coloneqq\\{i\in\\{1,\ldots,I_{0}I_{N}\\}\nonscript\,|\allowbreak\nonscript\,\mathopen{}[\operatorname{vec}(\mathbf{S}_{\Omega})]_{i}=1\\}$.
Also, for any $I_{0}I_{N}\times 1$ vector $\mathbf{v}$, let
$\mathbf{v}|_{\operatorname{vec}\Omega}$ be the
$\lvert\operatorname{vec}\Omega\rvert\times 1$ sub-vector of $\mathbf{v}$
formed by the entries of $\mathbf{v}$ at indices $\operatorname{vec}\Omega$.

According to [108, Prop. 26.13], the desired solution satisfies
$\mathbf{X}_{*}\in\mathcal{A}_{\Omega}\Leftrightarrow\mathscr{S}_{\Omega}(\mathbf{Y})=\mathscr{S}_{\Omega}(\mathbf{X}_{*})$
and
$\nabla\mathcal{L}_{X}(\mathbf{X}_{*})\in(\mathcal{A}_{\Omega}-\mathscr{S}_{\Omega}(\mathbf{Y}))^{\perp}=\\{\mathbf{X}\in\mathbb{C}^{I_{0}\times
I_{N}}\nonscript\,|\allowbreak\nonscript\,\mathopen{}\mathbf{0}=\mathscr{S}_{\Omega}(\mathbf{X})\\}^{\perp}=\\{\mathbf{X}\in\mathbb{C}^{I_{0}\times
I_{N}}\nonscript\,|\allowbreak\nonscript\,\mathopen{}\mathbf{0}=\mathscr{S}_{\Omega^{\complement}}(\mathbf{X})\\}$,
where $\perp$ stands for the orthogonal complement of a linear subspace.
Hence,

$\displaystyle\nabla\mathcal{L}_{X}(\mathbf{X}_{*})$ $\displaystyle{}={}$
$\displaystyle(1+\tau_{X})\mathbf{X}_{*}+\lambda_{L}\mathbf{L}_{\epsilon}^{\beta}\mathbf{X}_{*}\mathbf{\Delta}\mathbf{\Delta}^{\intercal}-\hat{\mathbf{M}}^{(n)}-\tau_{X}\hat{\mathbf{X}}^{(n)}$
$\displaystyle=$
$\displaystyle(1+\tau_{X})\mathbf{X}_{*}+\lambda_{L}\mathbf{L}_{\epsilon}^{\beta}[\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})+\mathscr{S}_{\Omega}(\mathbf{X}_{*})]\mathbf{\Delta}\mathbf{\Delta}^{\intercal}$
$\displaystyle-\hat{\mathbf{M}}^{(n)}-\tau_{X}\hat{\mathbf{X}}^{(n)}$
$\displaystyle=$
$\displaystyle(1+\tau_{X})\mathbf{X}_{*}+\lambda_{L}\mathbf{L}_{\epsilon}^{\beta}\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})\mathbf{\Delta}\mathbf{\Delta}^{\intercal}$
$\displaystyle+\lambda_{L}\mathbf{L}_{\epsilon}^{\beta}\mathscr{S}_{\Omega}(\mathbf{Y})\mathbf{\Delta}\mathbf{\Delta}^{\intercal}-\hat{\mathbf{M}}^{(n)}-\tau_{X}\hat{\mathbf{X}}^{(n)}\,.$

The requirement
$\mathscr{S}_{\Omega^{\complement}}(\nabla\mathcal{L}_{X}(\mathbf{X}_{*}))=\mathbf{0}$
yields

$\displaystyle(1+\tau_{X})\mathbf{S}_{\Omega^{\complement}}\odot\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})+\lambda_{L}\mathbf{S}_{\Omega^{\complement}}\odot\mathbf{L}_{\epsilon}^{\beta}\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})\mathbf{\Delta}\mathbf{\Delta}^{\intercal}$
$\displaystyle=-\lambda_{L}\mathbf{S}_{\Omega^{\complement}}\odot\mathbf{L}_{\epsilon}^{\beta}\mathscr{S}_{\Omega}(\mathbf{Y})\mathbf{\Delta}\mathbf{\Delta}^{\intercal}+\mathbf{S}_{\Omega^{\complement}}\odot(\hat{\mathbf{M}}^{(n)}+\tau_{X}\hat{\mathbf{X}}^{(n)})\,,$

and after vectorizing by $\operatorname{vec}(\cdot)$,

$\displaystyle(1+\tau_{X})\operatorname{diag}(\operatorname{vec}\mathbf{S}_{\Omega^{\complement}})\cdot\operatorname{vec}\left(\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})\right)$
$\displaystyle+\lambda_{L}\operatorname{diag}(\operatorname{vec}\mathbf{S}_{\Omega^{\complement}})\cdot(\mathbf{\Delta}\mathbf{\Delta}^{\intercal}\otimes\mathbf{L}_{\epsilon}^{\beta})\operatorname{vec}\left(\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})\right)$
$\displaystyle{}={}-\lambda_{L}\operatorname{diag}(\operatorname{vec}\mathbf{S}_{\Omega^{\complement}})\operatorname{vec}(\mathbf{L}_{\epsilon}^{\beta}\mathscr{S}_{\Omega}(\mathbf{Y})\mathbf{\Delta}\mathbf{\Delta}^{\intercal})$
$\displaystyle\hphantom{{}={}}+\operatorname{diag}(\operatorname{vec}\mathbf{S}_{\Omega^{\complement}})\operatorname{vec}(\hat{\mathbf{M}}^{(n)}+\tau_{X}\hat{\mathbf{X}}^{(n)})\,,$

where $\otimes$ denotes the Kronecker product. Consequently, by keeping only
the entries at indices $\operatorname{vec}(\Omega^{\complement})$,

$\displaystyle(1+\tau_{X})\operatorname{vec}\left(\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})\right)|_{\operatorname{vec}(\Omega^{\complement})}$
$\displaystyle+\lambda_{L}(\mathbf{\Delta}\mathbf{\Delta}^{\intercal}\otimes\mathbf{L}_{\epsilon}^{\beta})|_{\operatorname{vec}(\Omega^{\complement})\times\operatorname{vec}(\Omega^{\complement})}\cdot\operatorname{vec}\left(\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})\right)|_{\operatorname{vec}(\Omega^{\complement})}$
$\displaystyle{}={}-\lambda_{L}\operatorname{vec}(\mathbf{L}_{\epsilon}^{\beta}\mathscr{S}_{\Omega}(\mathbf{Y})\mathbf{\Delta}\mathbf{\Delta}^{\intercal})|_{\operatorname{vec}(\Omega^{\complement})}$
$\displaystyle\hphantom{{}={}}+\operatorname{vec}(\hat{\mathbf{M}}^{(n)}+\tau_{X}\hat{\mathbf{X}}^{(n)})|_{\operatorname{vec}(\Omega^{\complement})}\,,$

where
$(\mathbf{\Delta}\mathbf{\Delta}^{\intercal}\otimes\mathbf{L}_{\epsilon}^{\beta})|_{\operatorname{vec}(\Omega^{\complement})\times\operatorname{vec}(\Omega^{\complement})}$
is the submatrix extracted from
$\mathbf{\Delta}\mathbf{\Delta}^{\intercal}\otimes\mathbf{L}_{\epsilon}^{\beta}$
by keeping only the rows and columns indicated by
$\operatorname{vec}(\Omega^{\complement})$. It can be readily verified that

$\displaystyle\operatorname{vec}\bigl{(}\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})\bigr{)}|_{\operatorname{vec}(\Omega^{\complement})}$
$\displaystyle{}={}$ $\displaystyle\Bigl{(}(1$
$\displaystyle+\tau_{X})\mathbf{I}_{|\operatorname{vec}(\Omega^{\complement})\rvert}+\lambda_{L}(\mathbf{\Delta}\mathbf{\Delta}^{\intercal}\otimes\mathbf{L}_{\epsilon}^{\beta})|_{\operatorname{vec}(\Omega^{\complement})\times\operatorname{vec}(\Omega^{\complement})}\Bigr{)}^{-1}$
$\displaystyle\cdot\Bigl{(}$
$\displaystyle-\lambda_{L}\operatorname{vec}(\mathbf{L}_{\epsilon}^{\beta}\mathscr{S}_{\Omega}(\mathbf{Y})\mathbf{\Delta}\mathbf{\Delta}^{\intercal})|_{\operatorname{vec}(\Omega^{\complement})}$
$\displaystyle+\operatorname{vec}(\hat{\mathbf{M}}^{(n)}+\tau_{X}\hat{\mathbf{X}}^{(n)})|_{\operatorname{vec}(\Omega^{\complement})}\Bigr{)}\,.$
(19)

Recovering $\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})$ from (19) is a
straightforward task: take an all-zero $I_{0}I_{N}\times 1$ vector, replace
its zero entries at positions $\operatorname{vec}(\Omega^{\complement})$ by
the corresponding entries of (19), and then apply
$\operatorname{vec}^{-1}(\cdot)$ to recover
$\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})$. Finally,

$\displaystyle\hat{\mathbf{X}}^{(n+1/2)}=\mathbf{X}_{*}$
$\displaystyle=\mathscr{S}_{\Omega}(\mathbf{X}_{*})+\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})$
$\displaystyle=\mathscr{S}_{\Omega}(\mathbf{Y})+\mathscr{S}_{\Omega^{\complement}}(\mathbf{X}_{*})\,.$

###  C-B Solving (14b)

For compact notations, let
$\hat{\mathbfcal{D}}_{Q+1}^{(n)}\coloneqq\mathbfcal{K}$ and
$\hat{\mathbfcal{D}}_{Q+2}^{(n)}\coloneqq\hat{\mathbfcal{B}}^{(n)}$. For any
$q\in\\{2,\ldots,Q\\}$, define

$\displaystyle\hat{\mathbf{L}}_{q}^{(n)}$
$\displaystyle\coloneqq\prod\nolimits_{q^{\prime}=1}^{q-1}\hat{\mathbfcal{D}}_{q^{\prime}}^{(n)}\,,$
$\displaystyle\hat{\mathbf{R}}_{q}^{(n)}$
$\displaystyle\coloneqq\prod\nolimits_{q^{\prime}=q+1}^{Q+2}\hat{\mathbfcal{D}}_{q^{\prime}}^{(n)}\,,$

while $\hat{\mathbf{L}}_{1}^{(n)}\coloneqq\mathbf{I}_{I_{0}}$ and
$\hat{\mathbf{R}}_{1}^{(n)}\coloneqq\prod_{q^{\prime}=2}^{Q+2}\hat{\mathbfcal{D}}_{q^{\prime}}^{(n)}$,
so that the gradient of the loss in (14b) takes the form

$\displaystyle\nabla\mathcal{L}_{D}(\mathbfcal{D}_{q})$ $\displaystyle{}={}$
$\displaystyle\hat{\mathbf{L}}_{q}^{(n)}\mathbfcal{D}_{q}\hat{\mathbf{R}}_{q}^{(n)}-\hat{\mathbf{X}}^{(n)}$
$\displaystyle+\lambda_{2}\mathbfcal{D}_{q}+\tau_{D}(\mathbfcal{D}_{q}-\hat{\mathbfcal{D}}_{q}^{(n)})\,.$

The solution $\hat{\mathbfcal{D}}_{q}^{(n+1/2)}$ to (14b) takes a closed form,
and the way to obtain it follows the steps of Section C-A, where the affine
constraint to amount for the block-diagonal form of $\mathbfcal{D}_{q}$, for
$q\in\\{2,\ldots,Q\\}$, becomes
$\mathcal{A}_{\Omega}\coloneqq\\{\mathbfcal{D}_{q}\in\mathbb{C}^{d_{q-1}M\times
d_{q}M}{\mathbin{}\mid\mathbin{}}\mathbf{0}=\mathscr{S}_{\Omega}(\mathbfcal{D}_{q})\\}$,
with the index set $\Omega$ comprising all pairs of indices positioned off the
diagonal blocks. In the case of $\mathbfcal{D}_{1}$, no affine constraint (no
block-diagonal structure) is necessary to solve (14b). Solving
$\nabla\mathcal{L}_{D}(\mathbfcal{D}_{1})=\mathbf{0}$ for $\mathbfcal{D}_{1}$
will do the job.

##  Appendix D Solving the dMRI inverse problem

This section expands the details of the dMRI inverse problem (15). Similarly
to the discussion in Section IV-B, and upon defining the following tuple of
estimates $\forall n\in\mathbb{N},\forall k\in\\{0,1\\}$,

$\displaystyle\hat{\mathbfcal{O}}^{(n+k/2)}\coloneqq({}$
$\displaystyle{}\hat{\mathbf{X}}^{(n+k/2)},\hat{\mathbf{Z}}^{(n+k/2)},\hat{\mathbfcal{D}}_{1}^{(n+k/2)},$
$\displaystyle\ldots,\hat{\mathbfcal{D}}_{Q}^{(n+k/2)},\hat{\mathbfcal{B}}^{(n+k/2)})\,,$
(20)

Algorithm 1 solves (15) via computing solutions of the following convex sub-
tasks:

$\displaystyle\hat{\mathbf{X}}^{(n+1/2)}$
$\displaystyle\in\arg\min_{{\mathbf{X}}}{}$ $\displaystyle{}\
\tfrac{1}{2}\lVert{\mathbf{X}}-\hat{\mathbfcal{D}}_{1}^{(n)}\cdots\hat{\mathbfcal{D}}_{Q}^{(n)}\mathbfcal{K}\hat{\mathbfcal{B}}^{(n)}\rVert_{\textnormal{F}}^{2}$
$\displaystyle{}\
+\tfrac{\lambda_{2}}{2}\lVert\hat{\mathbf{Z}}^{(n)}-\mathscr{F}_{\text{t}}({\mathbf{X}})\rVert_{\textnormal{F}}^{2}$
$\displaystyle{}\
+\tfrac{\tau_{X}}{2}\lVert\mathbf{X}-\hat{\mathbf{X}}^{(n)}\rVert_{\textnormal{F}}^{2}$
s.to $\displaystyle{}\
\mathscr{S}_{\Omega}(\mathbf{Y})=\mathscr{S}_{\Omega}\mathscr{F}({\mathbf{X}})\,,$
(21a) $\displaystyle\hat{\mathbf{Z}}^{(n+1/2)}$
$\displaystyle\in\arg\min_{\mathbf{Z}}{}$ $\displaystyle{}\
\tfrac{\lambda_{2}}{2}\lVert\mathbf{Z}-\mathscr{F}_{\text{t}}(\hat{\mathbf{X}}^{(n)})\rVert^{2}_{\textnormal{F}}+\lambda_{3}\lVert\mathbf{Z}\rVert_{1}$
$\displaystyle{}\
+\tfrac{\tau_{Z}}{2}\lVert\mathbf{Z}-\hat{\mathbf{Z}}^{(n)}\rVert_{\textnormal{F}}^{2}\,,$
(21b) $\displaystyle\hat{\mathbfcal{D}}_{q}^{(n+1/2)}{}$
$\displaystyle{}\in\arg\min_{\mathbfcal{D}_{q}}{}$ $\displaystyle{}\
\tfrac{1}{2}\lVert\hat{\mathbf{X}}^{(n)}-\hat{\mathbfcal{D}}_{1}^{(n)}\cdots\mathbfcal{D}_{q}\cdots\hat{\mathbfcal{D}}_{Q}^{(n)}\mathbfcal{K}\hat{\mathbfcal{B}}^{(n)}\rVert_{\textnormal{F}}^{2}$
$\displaystyle{}\
+\tfrac{\lambda_{4}}{2}\lVert\mathbfcal{D}_{q}\rVert_{\textnormal{F}}^{2}+\tfrac{\tau_{D}}{2}\lVert\mathbfcal{D}_{q}-\hat{\mathbfcal{D}}_{q}^{(n)}\rVert_{\textnormal{F}}^{2}\,,$
s.to $\displaystyle{}\ \mathbfcal{D}_{q}\ \textnormal{is block diagonal}\
\forall q\in\\{2,\ldots,Q\\}\,,$ (21c)
$\displaystyle\hat{\mathbfcal{B}}^{(n+1/2)}$
$\displaystyle\in\arg\min_{\mathbfcal{B}}{}$ $\displaystyle{}\
\tfrac{1}{2}\lVert\hat{\mathbf{X}}^{(n)}-\hat{\mathbfcal{D}}_{1}^{(n)}\cdots\hat{\mathbfcal{D}}_{Q}^{(n)}\mathbfcal{K}\mathbfcal{B}\rVert_{\textnormal{F}}^{2}$
$\displaystyle{}\
+\lambda_{1}\lVert\mathbfcal{B}\rVert_{1}+\tfrac{\tau_{B}}{2}\lVert\mathbfcal{B}-\hat{\mathbfcal{B}}^{(n)}\rVert_{\textnormal{F}}^{2}$
s.to $\displaystyle{}\
\mathbf{1}_{N_{\mathit{l}}}^{{\mathsf{H}}}\mathbf{B}_{m}=\mathbf{1}_{I_{N}}^{{\mathsf{H}}}\,,\forall
m\in\\{1,\ldots,M\\}\,.$ (21d)

Sub-task (21d) is composite convex with affine constraints and can be solved
by [79]; details are provided in [48]. Sub-task (21c) is similar to sub-task
(14b) of the TVGS case, so that its solution follows Section C-B. The
following are the solutions to sub-tasks (21a) and (21b).

###  D-A Solving (21a)

Following [48], the unique solution $\hat{\mathbf{X}}^{(n+1/2)}$ to (14a) is

$\displaystyle\hat{\mathbf{X}}^{(n+1/4)}$ $\displaystyle\coloneqq
C_{X}\left(\hat{\mathbf{M}}^{(n)}+\lambda_{2}I_{N}\mathscr{F}_{t}^{-1}(\hat{\mathbf{Z}}^{(n)})+\tau_{X}\hat{\mathbf{X}}^{(n)}\right)\,,$
$\displaystyle\hat{\mathbf{X}}^{(n+1/2)}$
$\displaystyle\coloneqq\mathscr{F}^{-1}\mathscr{S}_{\Omega}(\mathbf{Y})+\mathscr{F}^{-1}\mathscr{S}_{\Omega^{\complement}}\mathscr{F}(\hat{\mathbf{X}}^{(n+1/4)})\,,$

where $C_{X}\coloneqq 1/(1+\lambda_{2}I_{N}+\tau_{X})$.

###  D-B Solving (21b)

Consider a convex function $g(\cdot):\mathbb{C}^{m\times
n}\to\mathbb{R}\cup\\{+\infty\\}$ and $\lambda\in\mathbb{R}_{>0}$. The
proximal mapping $\operatorname{Prox}_{\lambda g}(\cdot):\mathbb{C}^{m\times
n}\to\mathbb{C}^{m\times n}$ is defined as $\operatorname{Prox}_{\lambda
g}(\mathbf{A})\coloneqq\arg\min_{\mathbf{B}\in\mathbb{C}^{m\times n}}\lambda
g(\mathbf{B})+(1/2)\lVert\mathbf{A}-\mathbf{B}\rVert_{\text{F}}^{2}$ [108].
When $g$ is the $\ell_{1}$-norm $\lVert\cdot\rVert_{1}$, the $(i,j)$th entry
of $\operatorname{Prox}_{\lambda\lVert\cdot\rVert_{1}}(\mathbf{A})$ is
obtained by the soft-thresholding rule [108]:

$\displaystyle[\operatorname{Prox}_{\lambda\lVert\cdot\rVert_{1}}(\mathbf{A})]_{ij}=[\mathbf{A}]_{ij}\left(1-\tfrac{\lambda}{\max\\{\lambda,\lvert[\mathbf{A}]_{ij}\rvert\\}}\right)\,.$
(22)

Instead of $\operatorname{Prox}_{\lambda\lVert\cdot\rVert_{1}}(\mathbf{A})$,
it is more common to use notation $\textnormal{Soft}(\mathbf{A},\lambda)$. As
such, the solution $\hat{\mathbf{Z}}^{(n+1/2)}$ of (21b) is provided by

$\displaystyle\hat{\mathbf{Z}}^{(n+1/2)}=\text{Soft}[\mathscr{F}_{t}(\hat{\mathbf{X}}^{(n)})+\tfrac{\tau_{Z}}{\lambda_{2}}\hat{\mathbf{Z}}^{(n)},\tfrac{\lambda_{3}}{\lambda_{2}}]\,.$

##  Appendix E More results for TVGS recovery

This section provides additional figures and tables containing experimental
results of TVGS recovery (cf., Section IV-C).

TABLE V: Categorization of TVGS recovery methods. Method | Factorization/structured | Low-rank | Kernel-based
---|---|---|---
GraphTRSS [77] | ✗ | ✗ | ✗
KRG [55] | ✓ | ✗ | ✓
KGL [54] | ✓ | ✗ | ✓
LRDS [70] | ✗ | ✓ | ✗
MMF [14] | ✓ | ✓ | ✗
NBP [30] | ✓ | ✓ | ✓
MultiL-KRIM | ✓ | ✓ | ✓
Figure 6: Sensitivity of MAE to variation of parameters $d_{1}$ and kernel
choice in the context of (D3, P1, Nav1, L1) at 10% sampling ratio.

###  E-A Performance metrics

Recall that the original TVGS data is $\mathbf{Y}$ and the reconstructed data
is $\mathbf{X}$, which are both matrices of size $I_{0}\times I_{N}$. The
error metrics are defined as follows,

RMSE
$\displaystyle\coloneqq\frac{\lVert\mathbf{X}-\mathbf{Y}\rVert_{\textnormal{F}}}{\sqrt{I_{0}\cdot
I_{N}}}\,,$ MAE
$\displaystyle\coloneqq\frac{\lVert\mathbf{X}-\mathbf{Y}\rVert_{1}}{I_{0}\cdot
I_{N}}\,,$ MAPE $\displaystyle\coloneqq\frac{1}{I_{0}\cdot
I_{N}}\sum\nolimits_{i=1}^{I_{0}}\sum\nolimits_{j=1}^{I_{N}}\left\lvert\frac{x_{ij}-y_{ij}}{y_{ij}}\right\rvert\,.$

###  E-B Navigator- and landmark-point analysis

To examine the effect of different navigator-data formation and landmark-point
selection strategies, tests are conducted on the (D1, P1,$\cdot$,$\cdot$)
scenario. Table VI shows the error metrics of MultiL-KRIM under scenario (D1,
P1,$\cdot$, L1) when using different navigator-data formations. Averaged over
all sampling ratios, navigator data created from snapshots (Nav1) have
advantage over other formations, while those formed by nodes signals (Nav2)
produce highest errors. On the other hand, Table VII depicts variants of
landmark-point selection strategies under scenario (D1, P1, Nav1,$\cdot$).
Overall, there is no significant difference when using maxmin or clustering
methods for landmark-point selection.

TABLE VI: Performance on different navigators formation strategies under
scenario (D1, P1, _, L1).

Sampling ratio | Metric | Nav1 | Nav2 | Nav3 | Nav4
---|---|---|---|---|---
| MAE | 0.3004 | 0.3562 | 0.3350 | 0.3130
0.1 | RMSE | 0.4211 | 0.5065 | 0.5055 | 0.4219
| MAPE | 0.0296 | 0.0299 | 0.0283 | 0.0271
| MAE | 0.1732 | 0.1794 | 0.1883 | 0.1811
0.2 | RMSE | 0.2741 | 0.2899 | 0.3046 | 0.2872
| MAPE | 0.0175 | 0.0195 | 0.0200 | 0.0194
| MAE | 0.1176 | 0.1313 | 0.1213 | 0.1303
0.3 | RMSE | 0.1861 | 0.2074 | 0.1916 | 0.2045
| MAPE | 0.0137 | 0.0157 | 0.0150 | 0.0158
| MAE | 0.0824 | 0.0865 | 0.0834 | 0.0812
0.4 | RMSE | 0.1607 | 0.1699 | 0.1643 | 0.1635
| MAPE | 0.0082 | 0.0089 | 0.0082 | 0.0084
| MAE | 0.0562 | 0.0638 | 0.0585 | 0.0608
0.5 | RMSE | 0.1155 | 0.0994 | 0.0998 | 0.0949
| MAPE | 0.0050 | 0.0059 | 0.0050 | 0.0058
| MAE | 0.1460 | 0.1635 | 0.1573 | 0.1533
Average | RMSE | 0.2315 | 0.2546 | 0.2532 | 0.2344
| MAPE | 0.0148 | 0.0160 | 0.0153 | 0.0153

TABLE VII: Performance on different landmark points selection strategies under
scenario (D1, P1, Nav1, _).

Sampling ratio | Metric | L1 | L2 | L3
---|---|---|---|---
| MAE | 0.3004 | 0.3134 | 0.3290
0.1 | RMSE | 0.4211 | 0.4180 | 0.4587
| MAPE | 0.0296 | 0.0285 | 0.0291
| MAE | 0.1732 | 0.1785 | 0.1796
0.2 | RMSE | 0.2741 | 0.2847 | 0.2954
| MAPE | 0.0175 | 0.0179 | 0.0190
| MAE | 0.1176 | 0.1193 | 0.1166
0.3 | RMSE | 0.1861 | 0.1885 | 0.1848
| MAPE | 0.0137 | 0.0144 | 0.0144
| MAE | 0.0824 | 0.0819 | 0.0811
0.4 | RMSE | 0.1607 | 0.1605 | 0.1596
| MAPE | 0.0082 | 0.0083 | 0.0082
| MAE | 0.0562 | 0.0565 | 0.0551
0.5 | RMSE | 0.1155 | 0.0998 | 0.0996
| MAPE | 0.0050 | 0.0048 | 0.0048
| MAE | 0.1460 | 0.1499 | 0.1523
Average | RMSE | 0.2315 | 0.2303 | 0.2396
| MAPE | 0.0148 | 0.0147 | 0.0151

TABLE VIII: Runtime comparison (in seconds) between LRDS and MultiL-KRIM in
dataset D3.

Sampling pattern | Sampling rate | LRDS | MultiL-KRIM | Speedup
---|---|---|---|---
P1 | 0.1 | 1,393 | 413 | 3.4x
0.2 | 1,468 | 490 | 3x
0.3 | 1,358 | 719 | 1.9x
0.4 | 1,356 | 330 | 4.1x
0.5 | 1,346 | 270 | 5x
P2 | 0.1 | 1,433 | 545 | 2.6x
0.2 | 1,430 | 465 | 3.1x
0.3 | 1,463 | 461 | 3.2x
0.4 | 1,487 | 428 | 3.5x
0.5 | 1,470 | 562 | 2.6x

TABLE IX: Number of parameter unknowns (9) $Q$ | $\\{d_{q}\\}_{q=1}^{Q-1}$ | $N_{Q}[M=1]$ | $N_{Q}[M=7]$
---|---|---|---
1 | _ | $11,677,680$ | $81,743,760$
2 | $\\{8\\}$ | $1,357,472$ | $9,502,304$
3 | $\\{2,8\\}$ | $358,704$ | $2,510,928$
TABLE X: Number of parameter unknowns of MultiL-KRIM and TDDIP Method | Number of parameter unknowns
---|---
MultiL-KRIM[$M=7,Q=2$] | $9,502,304$
MultiL-KRIM[$M=1,Q=2$] | $1,357,472$
MultiL-KRIM[$M=7,Q=3$] | $2,510,928$
MultiL-KRIM[$M=1,Q=3$] | $358,704$
TDDIP1 | $3,162,025$
TDDIP2 | $2,492,713$
TDDIP3 | $394,985$
TDDIP4 | $255,417$
TABLE XI: Performance of MultiL-KRIM and TDDIP on radial sampling
(acceleration rate: 16x)

Methods $\setminus$ Metrics | NRMSE | SSIM | HFEN | M1 | M2 | Time
---|---|---|---|---|---|---
tableMulti MultiL-KRIM[$M=7,Q=2$] | ${{0.0448}}$ | ${0.8680}$ | ${\mathbf{0.1023}}$ | $\mathbf{37.35}$ | $\mathbf{1.4\times 10^{6}}$ | 1.2hrs
tableSingle MultiL-KRIM[$M=1,Q=2$] | ${0.0465}$ | ${0.8618}$ | ${0.1305}$ | $37.26$ | $\mathbf{1.4\times 10^{6}}$ | 44min
tableMulti MultiL-KRIM[$M=7,Q=3$] | $0.0448$ | ${0.8675}$ | ${0.1030}$ | $\mathbf{37.35}$ | $\mathbf{1.4\times 10^{6}}$ | 1hr
tableSingle MultiL-KRIM[$M=1,Q=3$] | ${0.0465}$ | ${0.8618}$ | ${0.1305}$ | $37.26$ | $\mathbf{1.4\times 10^{6}}$ | 41min
TDDIP1 [27] | $\mathbf{0.0376}$ | $\mathbf{0.8896}$ | $0.1452$ | ${36.24}$ | $1.3\times 10^{6}$ | 3.9hrs
TDDIP2 [27] | ${0.0388}$ | ${0.8855}$ | $0.1464$ | ${35.14}$ | $1.3\times 10^{6}$ | 3.5hrs
TDDIP3 [27] | ${0.0415}$ | ${0.8691}$ | $0.1642$ | ${35.87}$ | $\mathbf{1.4\times 10^{6}}$ | 2hrs
TDDIP4 [27] | ${0.0422}$ | ${0.8680}$ | $0.1668$ | ${35.68}$ | $\mathbf{1.4\times 10^{6}}$ | 1.2hrs

Dataset | Sampling ratio | Metric | LRDS | GraphTRSS | KRG | KGL | MMF | NBP | MultiL-KRIM
---|---|---|---|---|---|---|---|---|---
Sea Temperature |  | MAE | 0.3450 | 0.5822 | 0.4487 | 0.4438 | 0.3338 | 0.3009 | 0.3004
0.1 | RMSE | 0.5615 | 0.8645 | 0.5906 | 0.5824 | 0.4773 | 0.4018 | 0.4211
| MAPE | 0.0327 | 0.0507 | 0.0430 | 0.0405 | 0.0277 | 0.0269 | 0.0296
|  | MAE | 0.1817 | 0.3200 | 0.2478 | 0.2352 | 0.1758 | 0.1726 | 0.1732
| 0.2 | RMSE | 0.2814 | 0.5168 | 0.4156 | 0.4049 | 0.2806 | 0.2756 | 0.2741
|  | MAPE | 0.0186 | 0.0295 | 0.0222 | 0.0211 | 0.0185 | 0.0182 | 0.0175
|  | MAE | 0.1185 | 0.2093 | 0.1669 | 0.1649 | 0.1199 | 0.1206 | 0.1176
| 0.3 | RMSE | 0.1867 | 0.3659 | 0.2880 | 0.2832 | 0.1901 | 0.1903 | 0.1861
|  | MAPE | 0.0131 | 0.0200 | 0.0171 | 0.0169 | 0.0140 | 0.0142 | 0.0137
|  | MAE | 0.0857 | 0.1400 | 0.1103 | 0.1153 | 0.0847 | 0.0840 | 0.0824
| 0.4 | RMSE | 0.1422 | 0.2584 | 0.2328 | 0.2438 | 0.1662 | 0.1643 | 0.1607
|  | MAPE | 0.0098 | 0.0145 | 0.0092 | 0.0093 | 0.0083 | 0.0079 | 0.0082
|  | MAE | 0.0621 | 0.0982 | 0.0749 | 0.0782 | 0.0602 | 0.0566 | 0.0562
| 0.5 | RMSE | 0.1079 | 0.1966 | 0.1236 | 0.1292 | 0.0941 | 0.0889 | 0.1155
|  | MAPE | 0.0066 | 0.0101 | 0.0061 | 0.0063 | 0.0052 | 0.0055 | 0.0050
PM2.5 Concentration |  | MAE | 2.586 | 2.7698 | 2.5111 | 2.5214 | 2.3684 | 2.2936 | 2.2759
0.1 | RMSE | 4.558 | 4.7418 | 4.8145 | 4.8593 | 4.3193 | 4.2110 | 4.1863
| MAPE | 0.4928 | 0.4558 | 0.4778 | 0.4815 | 0.3941 | 0.3758 | 0.4053
|  | MAE | 2.0248 | 2.1572 | 1.7916 | 1.8028 | 1.7688 | 1.7041 | 1.6893
| 0.2 | RMSE | 3.7212 | 3.6645 | 3.3176 | 3.3203 | 3.2665 | 3.2666 | 3.1932
|  | MAPE | 0.3413 | 0.3765 | 0.3003 | 0.2994 | 0.3028 | 0.2890 | 0.2898
|  | MAE | 1.6123 | 1.7004 | 1.5179 | 1.4923 | 1.4380 | 1.3917 | 1.3255
| 0.3 | RMSE | 3.10697 | 3.0457 | 3.1716 | 3.1646 | 3.0090 | 2.9376 | 2.6244
|  | MAPE | 0.2907 | 0.292 | 0.2440 | 0.2380 | 0.2336 | 0.2347 | 0.2354
|  | MAE | 1.2519 | 1.4574 | 1.1746 | 1.1553 | 1.1251 | 1.0994 | 1.0342
| 0.4 | RMSE | 2.5083 | 2.6089 | 2.5277 | 2.5188 | 2.4154 | 2.3277 | 2.2939
|  | MAPE | 0.2504 | 0.2554 | 0.1947 | 0.1912 | 0.1933 | 0.1911 | 0.1668
Sea Level Pressure |  | MAE | 432.3553 | 438.5799 | 436.2410 | 435.4914 | 409.7814 | 408.8739 | 393.5182
0.1 | RMSE | 653.6381 | 673.0891 | 675.2066 | 678.8780 | 619.8458 | 618.0385 | 577.7076
| MAPE | 0.0044 | 0.0043 | 0.0041 | 0.0041 | 0.0040 | 0.0040 | 0.0038
|  | MAE | 315.3191 | 340.9664 | 323.9920 | 321.1525 | 293.8725 | 290.9654 | 277.9168
| 0.2 | RMSE | 502.2517 | 558.7858 | 555.1179 | 523.1133 | 477.8935 | 481.9912 | 463.2686
|  | MAPE | 0.0031 | 0.0034 | 0.0032 | 0.0029 | 0.0028 | 0.0028 | 0.0027
|  | MAE | 254.9302 | 275.5929 | 264.8433 | 261.0184 | 220.3090 | 220.6122 | 208.6543
| 0.3 | RMSE | 418.5478 | 476.5155 | 443.2093 | 430.8483 | 377.7021 | 380.5985 | 363.3772
|  | MAPE | 0.0025 | 0.0027 | 0.0025 | 0.0023 | 0.0023 | 0.0023 | 0.0021
|  | MAE | 201.3499 | 224.1063 | 206.5154 | 206.9803 | 162.9623 | 163.7134 | 158.0415
| 0.4 | RMSE | 352.2565 | 412.6796 | 367.3505 | 355.3074 | 303.4321 | 311.2039 | 296.1789
|  | MAPE | 0.0021 | 0.0022 | 0.0019 | 0.0018 | 0.0016 | 0.0016 | 0.0016
|  | MAE | 154.1453 | 179.4160 | 169.0522 | 163.3182 | 130.1137 | 123.4560 | 119.3680
| 0.5 | RMSE | 296.6426 | 352.4861 | 319.3949 | 321.6328 | 290.1288 | 269.4014 | 278.5416
|  | MAPE | 0.0015 | 0.0018 | 0.0016 | 0.0015 | 0.0013 | 0.0012 | 0.0012

TABLE XII: Performance comparisons for Random Sampling (P1). Red-colored
numbers are lowest errors.

Dataset | Sampling ratio | Metric | LRDS | GraphTRSS | KRG | KGL | MMF | NBP | MultiL-KRIM
---|---|---|---|---|---|---|---|---|---
Sea Temperature |  | MAE | 1.1540 | 1.2054 | 1.1705 | 1.1617 | 1.0946 | 1.0850 | 0.9966
0.1 | RMSE | 1.6223 | 1.7232 | 1.6684 | 1.6596 | 1.4288 | 1.4154 | 1.3542
| MAPE | 0.1172 | 0.1247 | 0.1172 | 0.1171 | 0.1175 | 0.1164 | 0.1157
|  | MAE | 0.8621 | 0.9962 | 0.8935 | 0.8786 | 0.8161 | 0.7927 | 0.7762