---
abstract: 'Matrix completion and extrapolation (MCEX) are dealt with here over reproducing kernel Hilbert spaces (RKHSs) in order to account for prior information present in the available data. Aiming at a fast and low-complexity solver, the task is formulated as one of kernel ridge regression. The resultant MCEX algorithm can also afford online implementation, while the class of kernel functions also encompasses several existing approaches to MC with prior information. Numerical tests on synthetic and real datasets show that the novel approach is faster than widespread methods such as alternating least-squares (ALS) or stochastic gradient descent (SGD), and that the recovery error is reduced, especially when dealing with noisy data.'
author:
- |
    Pere Giménez-Febrer$^1$, Alba Pagès-Zamora$^1$, and Georgios B. Giannakis$^2$\
    $^1$SPCOM Group, Universitat Politècnica de Catalunya-Barcelona Tech, Spain\
    $^2$Dept. of ECE and Digital Technology Center, University of Minnesota, USA
title: 'Matrix completion and extrapolation via kernel regression [^1]'
---

Matrix completion, extrapolation, RKHS, kernel ridge regression, graphs, online learning

Introduction
============

With only a subset of its entries available, matrix completion (MC) amounts to recovering the unavailable entries by leveraging just the low-rank attribute of the matrix itself [@candes]. The relevant task arises in applications as diverse as image restoration [@ji2010robust], sensor networks [@yi2015partial], and recommender systems [@koren2009matrix]. To save power for instance, only a fraction of sensors may collect and transmit measurements to a fusion center, where the available spatio-temporal data can be organized in a matrix format, and the unavailable ones can be eventually interpolated via MC [@yi2015partial]. Similarly, collaborative filtering of ratings given by users to a small number of items are stored in a sparse matrix, and the objective is to predict their ratings for the rest of the items [@koren2009matrix].

Existing MC approaches rely on some form of rank minimization or low-rank matrix factorization. Specifically,  [@candes] proves that when MC is formulated as the minimization of the nuclear norm subject to the constraint that the observed entries remain unchanged, exact recovery is possible under mild assumptions; see also [@candes2010matrix] where reliable recovery from a few observations is established even in the presence of additive noise. Alternatively, [@koren2009matrix] replaces the nuclear norm by two low-rank factor matrices that are identified in order to recover the complete matrix.

While the low-rank assumption can be sufficient for reliable recovery, prior information about the unknown matrix can be also accounted to improve the completion outcome. Forms of prior information can include sparsity [@yi2015partial], local smoothness [@cheng2013stcdg], and interdependencies encoded by graphs [@kalofolias2014matrix; @chen; @rao2015collaborative; @ma2011recommender]. These approaches exploit the available similarity information or prior knowledge of the bases spanning the column or row spaces of the unknown matrix. In this regard, reproducing kernel Hilbert spaces (RKHSs) constitute a powerful tool for leveraging available prior information thanks to the kernel functions, which measure the similarity between pairs of points in an input space. Prompted by this, [@abernethy2006low; @bazerque2013; @zhou2012; @stock2018comparative] postulate that columns of the factor matrices belong to a pair of RKHSs spanned by their respective kernels. In doing so, a given structure or similarity between rows or columns is effected on the recovered matrix. Upon choosing a suitable kernel function, [@yi2015partial] as well as [@cheng2013stcdg; @kalofolias2014matrix; @chen; @rao2015collaborative; @ma2011recommender] can be cast into the RKHS framework. In addition to improving MC performance, kernel-based approaches also enable extrapolation of rows and columns, even when all their entries are missing - a task impossible by the standard MC approaches in e.g. [@candes] and [@koren2009matrix].

One major hurdle in MC is the computational cost as the matrix size grows. In its formulation as a rank minimization task, MC can be solved via semidefinite programming [@candes], or proximal gradient minimization [@ma2011fixed; @cai2010singular; @chen; @gimenez], which entails a singular value decomposition of the recovered matrix per iteration. Instead, algorithms with lower computational cost are available for the bi-convex formulation based on matrix factorization [@koren2009matrix]. These commonly rely on iterative minimization schemes such as alternating least-squares (ALS) [@hastie2015matrix; @jain2013low] or stochastic gradient descent (SGD) [@gemulla2011large; @zhou2012]. With regard to kernel-based MC, the corresponding algorithms rely on alternating convex minimization and semidefinite programming [@abernethy2006low], block coordinate descent [@bazerque2013], and SGD [@zhou2012]. However, algorithms based on alternating minimization only converge to the minimum after infinite iterations. In addition, existing kernel-based algorithms adopt a specific sampling pattern or do not effectively make use of the Representer Theorem for RKHSs that will turn out to be valuable in further reducing the complexity, especially when the number of observed entries is small.

The present contribution offers an RKHS-based approach to MCEX that also unifies and broadens the scope of MC approaches, while offering reduced complexity algorithms that scale well with the data size. Specifically, we develop a novel MC solver via kernel ridge regression as a convex alternative to the nonconvex factorization-based formulation that offers a closed-form solution. Through an explicit sampling matrix, the proposed method offers an encompassing sampling pattern, which further enables the derivation of upper bounds on the mean-square error. Moreover, an approximate solution to our MCEX regression formulation is developed that also enables online implementation using SGD. Finally, means of incorporating prior information through kernels is discussed in the RKHS framework.

The rest of the paper paper is organized as follows. Section II outlines the RKHS formulation and the kernel regression task. Section III unifies the existing methods for MC under the RKHS umbrella, while Section IV introduces our proposed Kronecker kernel MCEX (KKMCEX) approach. Section V develops our ridge regression MCEX (RRMCEX) algorithm, an accelerated version of KKMCEX, and its online variant. Section VI deals with the construction of kernel matrices. Finally, Section VII presents numerical tests, and Section VIII concludes the paper.

**Notation.** Boldface lower case fonts denote column vectors, and boldface uppercase fonts denote matrices. The $(i,j)$th entry of matrix $\A$ is $\A_{i,j}$, and the $i^{\text{th}}$ entry of vector $\bm a$ is $\bm a_i$. Superscripts $^T$ and $^\dagger$ denote transpose and pseudoinverse, respectively; while hat $\,\hat{}\,$ is used for estimates. Matrix $\F\in\mathcal{H}$ means that its columns belong to a vector space $\mathcal{H}$. The symbols $\I$ and $\bm 1$ stand for the identity matrix and the all-ones vector of appropriate size, specified by the context. The trace operator is $\text{Tr}(\cdot)$, the function eig($\A$) returns the diagonal eigenvalue matrix of $\A$ ordered in ascending order, and $\lambda_k(\A)$ denotes the $k^{\text{th}}$ eigenvalue of $\A$ with $\lambda_k(\A)\leq\lambda_{k+1}(\A)$.

Preliminaries {#sec:background}
=============

Consider a set of $N$ input-measurement pairs $\{(x_i,m_i)\}^N_{i=1}$ in $\mathcal{X}\times\mathbb{R}$, where $\mathcal{X}:=\{x_1,\ldots,x_N\}$ is the input space, $\mathbb{R}$ denotes the set of real numbers, and measurements obey the model $$\label{eq:sigmodel}
  m_i = f(x_i) + e_i$$ where $f:\mathcal{X}\rightarrow\mathbb{R}$ is an unknown function and $e_i\in\mathbb{R}$ is noise. We assume this function belongs to an RKHS $$\label{eq:hilb}
  \mathcal{H}_x:=\{f:f(x_i) = \sum^N_{j=1} \alpha_j \kapx(x_i,x_j), \:\:\: \alpha_j \in \mathbb{R}\}$$ where $\kapx:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ is the kernel function that spans $\mathcal{H}_x$, and $\{\alpha_i\}^N_{i=1}$ are weight coefficients. An RKHS is a complete linear space endowed with an inner product that satisfies the reproducing property [@shawe]. If $\langle f,f'\rangle_{\mathcal{H}_x}$ denotes the inner product in $\mathcal{H}_x$ between functions $f$ and $f'$, the reproducing property states that $f(x_i) = \langle f,\kapx(\cdot,x_i)\rangle_{\mathcal{H}_x}$; that is, $f$ in $\mathcal{H}_x$ can be evaluated at $x_i$ by taking the inner product between $f$ and $\kapx(\cdot,x_i)$. With $\{\alpha_i\}^N_{i=1}$ and $\{\alpha'_i\}^N_{i=1}$ denoting the coefficients of $f$ and $f'$ in (\[eq:hilb\]) respectively, we have $\langle f,f'\rangle_{\mathcal{H}_x} := \sum_{i=1}^N\sum_{j=1}^N \alpha_i\alpha_j' \kapx(x_i,x_j)$; that is, $$\begin{aligned}
  \langle f,f'\rangle_{\mathcal{H}_x} = \balp^T\K_x\balp'\label{eq:innerprod}\end{aligned}$$ where $\balp := [\alpha_1,\ldots,\alpha_N]^T$, $\balp' := [\alpha_1',\ldots,\alpha_N']^T$ and $(\K_x)_{i,j}:=\kapx(x_i,x_j)$. In order for $\langle \cdot,\cdot \rangle_{\mathcal{H}_x}$ in (\[eq:innerprod\]) to be an inner product, $\kapx$ must be symmetric and semipositive definite, meaning $\langle f,f\rangle_{\mathcal{H}_x} \geq 0\: \forall f\in\mathcal{H}_x$. As a consequence, $\K_x$ will be symmetric positive semidefinite since $\balp^T\K_x\balp \geq 0 \:\: \forall\, \balp \in \mathbb{R}^N$.

While $\kapx$ is usually interpreted as a measure of similarity between two elements in $\mathcal{X}$, it can also be seen as the inner product of corresponding two elements in feature space $\mathcal{F}$ to which $\mathcal{X}$ can be mapped using function $\phi_x:\mathcal{X}\rightarrow\mathcal{F}$. Formally, we write $$\label{eq:kerdot}
    \kapx(x_i,x_j)=\langle\phi_x(x_i),\phi_x(x_j)\rangle_{\mathcal{F}}.$$ Function $\phi_x$ is referred to as feature map, and its choice depends on the application. For an input space of text files, for example, the files could be mapped to a feature vector that tracks the number of words, lines, and blank spaces in the file. Since $\phi_x$ can potentially have infinite dimension, evaluating the kernel using  might be prohibitively expensive. This motivates specifying the kernel through a similarity function in $\mathcal{X}$, which bypasses the explicit computation of the inner product in $\mathcal{F}$. Typical examples include the Gaussian kernel $\kapx(x_i,x_j) = \text{exp}\{-{\left|\left|x_i-x_j\right|\right|^2_2}/(2\eta)\}$ with $\eta$ being a free parameter, and the polynomial kernel [@friedman2001elements]. In certain cases however, it is difficult to obtain the kernel similarity function on the input space. Such cases include metric input spaces with misses (as in MC), and non-metric spaces. The alternative to both is deriving the kernel from prior information. For instance, if we have a graph connecting the points in $\mathcal{X}$, a kernel can be obtained from the graph Laplacian [@romero].

Having introduced the basics of RKHS, we proceed with the kernel regression task, where given $\{m_i\}^N_{i=1}$ we seek to obtain $$\label{eq:rrf}
    \hat{f} = \operatorname*{\arg\,\min}_{f\in\mathcal{H}_x} {1\over N}\sum_{i=1}^N l(m_i,f(x_i)) +  \mu' ||f||^2_{\mathcal{H}_x}$$ with $l(\cdot)$ denoting the loss, $\mu' \in \mathbb{R}^+$ the regularization parameter, and $||f||_{\mathcal{H}_x}:=\langle f,f\rangle_{\mathcal{H}_x}$ the norm induced by the inner product in (\[eq:innerprod\]). We will henceforth focus on the square loss $l(m_i,f(x_i)) := (m_i - f(x_i))^2$. Using $\K_x$, consider without loss of generality expressing the vector $\f:=[f(x_1),\ldots,f(x_N)]^T$ as $\f=\K_x\balp$, where $\balp:=[\alpha_1,\ldots,\alpha_N]^T$. Using the latter in the square loss, (\[eq:rrf\]) boils down to a kernel ridge regression (KRR) problem that can be solved by estimating $\balp$ as $$\label{eq:rralp}
  \hat{\balp} = \operatorname*{\arg\,\min}_{\balp\in\mathbb{R}^N} {\left|\left|\m - \K_x \balp\right|\right|^2_2} +\mu \mspace{2mu} \balp^T \K_x \balp$$ where $\m := [m_1,\ldots,m_N]^T$ and $\mu=N\mu'$. The weights can be found in closed form as $$\label{eq:rralpest}
\hat{\balp} = (\K_x+\mu\I)^{-1}\m$$ and the estimate of the sought function is obtained as $\hat{\f}=\K_x\hat{\balp}$.

Kernel-based MCEX {#sec:mc}
=================

Matrix completion considers $\F\in\mathbb{R}^{N\times L}$ of rank $r$ observed through a $N\times L$ matrix of noisy observations $$\label{eq:matm}
    \M = P_{\Omega}(\F+\bm \E)$$ where $\Omega\subseteq\{1,\ldots,N\}\times \{1,\ldots,L\}$ is the sampling set of cardinality $S=|\Omega|$ containing the indices of the observed entries; $P_{\Omega}(\cdot)$ is a projection operator that sets to zero the entries with index $(i,j)\notin \Omega$ and leaves the rest unchanged; and, $\bm\E\in\mathbb{R}^{N\times L}$ is a noise matrix. According to [@candes2010matrix], one can recover $\bm F$ from $\M$ with an error proportional to the magnitude of ${\left|\left|\E\right| \right|_\text{F}^2}$ by solving the following convex optimization problem: $$\begin{aligned}
\label{eq:nnormc}
    &\min_{\F\in{\mathbb{R}^{N\times L}}}
    & \mspace{-50mu}& \text{rank}(\F)  \nonumber \\
    & \text{subject to}
    & \mspace{-50mu}& {\left|\left|P_{\Omega}(\F - \M)\right| \right|_\text{F}^2} \leq \delta\end{aligned}$$ where ${\left|\left|\cdot\right| \right|_\text{F}}$ is the Frobenius norm, and we assume ${\left|\left|P_\Omega(\E)\right| \right|_\text{F}^2} \leq \delta$ for some $\delta > 0$. Since solving (\[eq:nnormc\]) is NP-hard, the nuclear norm ${\left|\left|\F\right|\right|_*} := \text{Tr}(\sqrt{\F^T\F})$ can be used to replace the rank to obtain the convex problem [@ma2011fixed; @chen] $$\label{eq:nuclmc}
  \min_{\F\in{\mathbb{R}^{N\times L}}} {\left|\left|P_\Omega(\M-\F)\right| \right|_\text{F}^2} + \mu{\left|\left|\F\right|\right|_*}.$$ Because $\F$ is low rank, it is always possible to factorize it as $\F=\W\H^T$, where $\W \in \mathbb{R}^{N\times p}$ and $\H \in \mathbb{R}^{L\times p}$ are the latent factor matrices with $p \geq r$. This factorization allows expressing the nuclear norm as [@srebro2005maximum] ${\left|\left|\F\right|\right|_*} = \min_{\F=\W\H^T} {1\over 2}\left({\left|\left|\W\right| \right|_\text{F}^2} + {\left|\left|\H\right| \right|_\text{F}^2}\right)$ which allows reformulating (\[eq:nuclmc\]) as $$\label{eq:mclag}
    \{\hat{\W}\!,\!\hat{\H}\} \!=\! \operatorname*{\arg\,\min}_{\substack{\W\in\mathbb{R}^{N\times p}\\\H\in\mathbb{R}^{L\times p}}} {\left|\left|P_\Omega(\!\M\!-\!\W\H^T)\right| \right|_\text{F}^2} \!+ \mu\!\left({\left|\left|\W\right| \right|_\text{F}^2} \!+\! {\left|\left|\H\right| \right|_\text{F}^2}\right)$$ and yields $\hat{\F}=\hat{\W}\hat{\H}^T$. While the solutions to (\[eq:nuclmc\]) and (\[eq:mclag\]) are equivalent when the rank of the matrix minimizing  is smaller than $p$[@hastie2015matrix], solving (\[eq:nuclmc\]) can be costlier since it involves the computation of the singular values of the matrix. On the other hand, since (\[eq:mclag\]) is bi-convex it can be solved by alternately optimizing $\W$ and $\H$, e.g. via ALS [@jain2013low] or SGD iterations [@gemulla2011large]. Moreover, leveraging the structure of (\[eq:mclag\]), it is also possible to optimize one row from each factor matrix at a time instead of updating the full factor matrices, which enables faster and also online and distributed implementations [@teflioudi2012].

Aiming at a kernel-based MCEX, we model the columns and rows of $\F$ as functions that belong to two different RKHSs. To this end, consider the input spaces $\mathcal{X}:=\{x_1,\ldots,x_N\}$ and $\mathcal{Y}:=\{y_1,\ldots,y_L\}$ for the column and row functions, respectively. In the user-movie ratings paradigm, $\mathcal{X}$ could be the set of users, and $\mathcal{Y}$ the set of movies. Then $\F:=[\f_1,\ldots,\f_L]$ is formed with columns $\f_l := [f_l(x_1),\ldots,f_l(x_N)]^T$ with $f_l:\mathcal{X}\rightarrow\mathbb{R}$. Likewise, we rewrite $\F := [\g_1,\ldots,\g_N]^T$, with rows $\g_n^T:= [g_n(y_1),\ldots,g_n(y_L)]$ and $g_n:\mathcal{Y}\rightarrow\mathbb{R}$ . We further assume that $f_l\in\mathcal{H}_x\: \forall l=1,\ldots,L$ and $g_n\in\mathcal{H}_y \: \forall n=1,\ldots,N$, where $$\label{eq:hilbrow}
  \mathcal{H}_x:=\{f:f(x_i) = \sum^N_{j=1} \alpha_j \kapx(x_i,x_j), \:\:\: \alpha_j \in \mathbb{R}\}$$ $$\label{eq:hilbcol}
  \mathcal{H}_y:=\{g:g(y_i) = \sum^L_{j=1} \beta_j \kapy(y_i,y_j), \:\:\: \beta_j \in \mathbb{R}\}$$ and $\kapx:\mathcal{X}\times\mathcal{X}\rightarrow \mathbb{R}$ and $\kapy:\mathcal{Y}\times\mathcal{Y}\rightarrow \mathbb{R}$ are the kernels forming $\K_x\in\mathbb{R}^{N\times N}$ and $\K_y\in\mathbb{R}^{L\times L}$, respectively. Since $\W$ and $\H$ span the column and row spaces of $\F$, their columns belong to $\mathcal{H}_x$ and $ \mathcal{H}_y$ as well. Thus, the $m^\text{th}$ column of $\W$ is $$\label{eq:vecw}
    \w_m:=[w_m(x_1),\ldots,w_m(x_N)]^T$$ where $w_m:\mathcal{X}\rightarrow\mathbb{R}$ and $w_m\in\mathcal{H}_x\:\forall m=1,\ldots,p$, and the $m^\text{th}$ column of $\H$ is $$\label{eq:vech}
\h_m:=[h_m(y_1),\ldots,h_m(y_L)]^T$$ where $h_m:\mathcal{Y}\rightarrow\mathbb{R}$ and $h_m\in\mathcal{H}_y\:\forall m=1,\ldots,p$. Hence, instead of simply promoting a small Frobenius norm for the factor matrices as in (\[eq:mclag\]), we can also promote smoothness on their respective RKHS. The kernel-based formulation in [@bazerque2013] estimates the factor matrices by solving $$\begin{aligned}
    \{\hat{\W},\hat{\H}\}=\operatorname*{\arg\,\min}_{\substack{\W\in\mathcal{H}_x\\\H\in\mathcal{H}_y}} \ &{\left|\left|P_\Omega(\M-\W\H^T)\right| \right|_\text{F}^2}  \label{eq:factker}\\& \mspace{-50mu}+ \mu{\text{Tr}(\W^T\K_x^{-1}\W)} + \mu{\text{Tr}(\H^T\K_y^{-1}\H)}.\nonumber\end{aligned}$$ Note that (\[eq:factker\]) is equivalent to (\[eq:mclag\]) for $\K_x=\I$ and $\K_y=\I$. Since the constraints $\W\in\mathcal{H}_x$ and $\H\in\mathcal{H}_y$ can be challenging to account for when solving (\[eq:factker\]), we can instead find the coefficients that generate $\W$ and $\H$ in their respective RKHSs in order to satisfy such constraints. Thus, if we expand $\W=\K_x\B$ and $\H=\K_y\C$, where $\B\in\mathbb{R}^{N\times p}$ and $\C\in\mathbb{R}^{L\times p}$ are coefficient matrices, (\[eq:factker\]) becomes $$\begin{aligned}
\label{eq:factkerw}
    \{\hat{\B},\hat{\C}\}=\operatorname*{\arg\,\min}_{\substack{\B\in\mathbb{R}^{N\times p}\\\C\in\mathbb{R}^{L\times p}}} & {\left|\left|P_\Omega(\M-\K_x\B\C^T\K_y)\right| \right|_\text{F}^2}  \\
    &+ \mu{\text{Tr}(\B^T\K_x\B)} + \mu{\text{Tr}(\C^T\K_y\C)}.\nonumber\end{aligned}$$ Nevertheless, with nonsingular kernel matrices, $\B$ and $\C$ can be found by solving and substituting $\hat{\B}=\K_x^{-1}\hat{\W}$ and $\hat{\C}=\K_y^{-1}\hat{\H}$ [@bazerque2013].

Alternating minimization schemes that solve the bi-linear MC formulation tends to the solution to the convex problem in the limit [@jain2013low], thus convergence to the global optimum is not guaranteed unless the number of iterations is infinite. Since algorithms for kernel-based MC [@bazerque2013] solving rely on such alternating minimization schemes, they lack convergence guarantees given finite iterations as well. In addition to that, their computational cost scales with the size of $\F$. On the other hand, online implementations have a lower cost [@zhou2012], but only guarantee convergence to a stationary point [@mardani2015]. In the ensuing section we develop a convex kernel-based reformulation of MCEX that enables a closed-form solver which purely exploits the extrapolation facilitated by the kernels. By casting aside the low-rank constraints, the computational complexity of our solver scales only with the number of observations while, according to our numerical tests, providing better performance. Moreover, we derive an online implementation that can be seamlessly extended to distributed operation.

Kronecker kernel MCEX {#sec:kmr}
=====================

In the previous section, we viewed the columns and rows of $\F$ as functions evaluated at the points of the input spaces $\mathcal{X}$ and $\mathcal{Y}$ in order to unify the state-of-the-art on MC using RKHSs. Instead, we now postulate entries of $\F$ as the output of a function lying on an RKHS evaluated at a tuple $(x_i,y_i) \in \mathcal{X} \times \mathcal{Y}$. Given the spaces $\mathcal{X}$ and $\mathcal{Y}$, consider the space $\mathcal{Z}:=\mathcal{X}\times\mathcal{Y}$ with cardinality $|\mathcal{Z}|=NL$ along with the two-dimensional function $v:\mathcal{Z}\rightarrow\mathbb{R}$ as $v(x_i,y_j)=f_j(x_i)$, which belongs to the RKHS $$\label{eq:hilbz}
  \mathcal{H}_z\!:=\!\{v\!:\!v(x_i,\!y_j) \!=\! \sum^N_{n=1}\!\sum^L_{l=1} \!\gamma_{n,l} \kapz((x_i,\!y_j),\!(x_n,\!y_l)), \: \:\; \gamma_{n,l} \!\in \!\mathbb{R}\}$$ with $\kapz:\mathcal{Z}\times\mathcal{Z}\rightarrow \mathbb{R}$. While one may choose any kernel to span $\mathcal{H}_z$, we will construct one adhering to our bilinear factorization $\F=\W\H^T$ whose $(i,j)^\text{th}$ entry yields $$\label{eq:zhw}
\F_{ij}=v(x_i,y_j) = \sum^p_{m=1}w_m(x_i)h_m(y_j)$$ with $w_m$ and $h_m$ functions capturing $m^{\text{th}}$ column vector of $\W$ and $\H$ as in (\[eq:vecw\]) and (\[eq:vech\]). Since $ w\in\mathcal{H}_x$ and $ h\in\mathcal{H}_y$, we can write $w_m(x_i) = \sum^N_{n=1} b_{n,m} \kapx(x_i,x_n)$ and $h_m(y_j) = \sum^L_{l=1} c_{l,m} \kapy(y_j,y_l)$, where $b_{n,m}$ and $c_{l,m}$ are the entries at $(n,m)$ and $(l,m)$ of the factor matrices $\B$ and $\C$ from , respectively. Therefore, (\[eq:zhw\]) can be rewritten as $$\begin{aligned}
v(x_i,y_j) &= \sum^p_{m=1} \sum^N_{n=1} b_{n,m}\kapx(x_i,x_n) \sum^L_{l=1} c_{l,m}\kapy(y_j,y_l) \nonumber \\
&= \sum^N_{n=1}\sum^L_{l=1} \left(\sum^p_{m=1} b_{n,m} c_{l,m} \right) \kapx(x_i,x_n)\kapy(y_j,y_l) \nonumber \\
&=  \sum^N_{n=1}\sum^L_{l=1} \gamma_{n,l} \kapz((x_i,y_j),(x_n,y_l))\label{eq:zfunc}\end{aligned}$$ where $\gamma_{n,l}=\sum^p_{m=1} b_{m,n} c_{m,l}$, and $\kapz((x_i,y_j),(x_n,y_l)) =\kapx(x_i,x_n)\kapy(y_j,y_l)$ since a product of kernels is itself a kernel [@friedman2001elements]. Using the latter,  can be written compactly as $$\label{eq:zkermat}
    v(x_i,y_j) = \k_{i,j}^T\gam$$ where $\gam:=[\gamma_{1,1},\gamma_{2,1},\ldots,\gamma_{N,1},\gamma_{1,2},\gamma_{2,2},\ldots,\gamma_{N,L}]^T$, and correspondingly, $$\begin{aligned}
 \label{eq:kykz}
    \k_{i,j}=&[\kapx(x_i,x_1)\kapy(y_j,y_1), \ldots, \kapx(x_i,x_N)\kapy(y_j,y_1),\nonumber \\ &\kapx(x_i,x_1)\kapy(y_j,y_2), \ldots, \kapx(x_i,x_N)\kapy(y_j,y_L)]^T \nonumber \\
    =& (\K_y)_{:,j}\otimes(\K_x)_{:,i}\end{aligned}$$ where a subscript $(:,j)$ denotes the $j^{\text{th}}$ column of a matrix, and we have used that $\K_x$ and $\K_y$ are symmetric matrices. In accordance with , the kernel matrix of $\mathcal{H}_z$ in  is $$\label{eq:Kz}
\K_z=\K_y\otimes\K_x.$$ Clearly, $\k_{i,j}$ in  can also be expressed as $\k_{i,j}=(\K_z)_{:,(j-1)N+i}$. This together with  implies that $$\begin{aligned}
  \v =[&v(x_1,y_1),v(x_2,y_1),\ldots,v(x_N,y_1),v(x_1,y_2),\\&v(x_2,y_2),\ldots,v(x_N,y_N)]^T \end{aligned}$$ can be expressed in matrix-vector form as $$\label{eq:zmod}
\v = \K_z\gam$$ or, equivalently, $\v=\text{vec}(\F)$. Note that entries of the kernel matrix are $(\K_z)_{i',j'} = \kapx(x_i,x_n)\kapy(y_j,y_l)$, where $n=j'\:\text{mod}\:N, i=i'\:\text{mod}\:N, l=\lceil {j'\over N}\rceil, \text{ and } j=\lceil {i'\over N}\rceil$.

Since the eigenvalues of $\K_z$ are the product of eigenvalues of $\K_y$ and $\K_x$, it follows that $\K_z$ is positive semidefinite and thus a valid kernel matrix. With the definition of the function $v$ and its vector form we have transformed the matrix of functions specifying $\F$ into a function that lies on the RKHS $\mathcal{H}_z$. Hence, we are ready to formulate MCEX as a kernel regression task for recovering $\v$ from the observed entries of $\m=\text{vec}(\M)$.

Given $ \{((x_i,y_j),m_{i,j})\}_{(i,j)\in\Omega}$ in $\mathcal{Z}\times\mathbb{R}$, our goal is to recover the function $v$ as $$\label{eq:lossr1}
  \hat{v}=\operatorname*{\arg\,\min}_{v\in\mathcal{H}_z} \sum_{(i,j)\in\Omega} (m_{i,j}-v(x_i,y_j))^2 + \mu ||v||^2_{\mathcal{H}_\mathcal{Z}}$$ where $||v||^2_{\mathcal{H}_\mathcal{Z}}:=\gam^T\K_z\gam$. Define next $\e := \text{vec}(\E)$ and $\mb=\S\m$, where $\S$ is an $S\times NL$ binary sampling matrix also used to specify the sampled noise vector $\bar{\e}=\S\e$. With these definitions and (\[eq:zmod\]), the model in (\[eq:matm\]) becomes $$\label{eq:obsvec}
  \mb = \S\v + \S\e = \S\K_z\gam + \bar{\e}$$ which can be solved to obtain $$\label{eq:mcker}
  \hgam = \operatorname*{\arg\,\min}_{\gam\in\mathbb{R}^{NL}} {\left|\left|\mb - \S\K_z\gam\right|\right|^2_2} + \mu\gam^T\K_z\gam$$ in closed form $$\label{eq:krr}
  \hgam = (\S^T\S\K_z+\mu\I)^{-1}\S^T\mb.$$ Since the size of $\K_z$ is $NL\times NL$, the inversion in (\[eq:krr\]) can be very computationally intensive. To alleviate this, we will leverage the Representer Theorem (see [@scholkopf2001generalized] for a formal proof), which allows us to reduce the number of degrees of freedom of the regression problem. In our setup, this theorem is as follows.

**Representer Theorem**. Given the set of input-observations pairs $\{(x_i,y_j),m_{i,j})\}_{(i,j)\in \Omega}$ in $\mathcal{Z}\times \mathbb{R}$ and the function $v$ as in (\[eq:zfunc\]), the solution to $$\label{eq:lossr}
  \operatorname*{\arg\,\min}_{v\in\mathcal{H}_z} \sum_{(i,j)\in\Omega} (m_{i,j}-v(x_i,y_j))^2 + \mu ||v||^2_{\mathcal{H}_\mathcal{Z}}$$ is an estimate $\hat{v}$ that satisfies $$\label{eq:zhat}
  \hat{v} = \sum_{(n,l)\in\Omega} {\tau}_{n,l} k_{z}((\cdot,\cdot),(x_n,y_l))$$ for some coefficients $\tau_{n,l}\in\mathbb{R},$ $\forall(n,l)\in\Omega$.

Theorem 1 asserts that $\hat{\gam}$ in (\[eq:mcker\]) satisfies $\hat{\gamma}_{n,l}=0 \; \forall \; (n,l) \notin \Omega$. Therefore, we only need to optimize $\{\gamma_{n,l} : (n,l)\in\Omega\}$ which correspond to the observed pairs. In fact, for our vector-based formulation, the Representer Theorem boils down to applying on  the matrix inversion lemma (MIL), which asserts the following.

\[lem:MIL\] **MIL** [@henderson1981]. Given matrices $\A,\U$ and $\V$ of conformal dimensions, with $\A$ invertible, it holds that $$\label{eq:mil}
   (\U\V+\A)^{-1}\U = \A^{-1}\U(\V\A^{-1}\U+\I).$$

With $\A=\mu\I, \U=\S^T$ and $\V=\K_z\S^T$, application of to   yields $$\label{eq:gamest}
  \hat{\gam} = \S^T(\S\K_z\S^T+\mu\I)^{-1}\mb.$$ Subsequently, we reconstruct $\v$ as $$\label{eq:zKKMCEX}
  \hat{\v}_K = \K_z\S^T(\S\K_z\S^T+\mu\I)^{-1}\mb$$ and we will henceforth refer to as the *Kronecker kernel MCEX* (KKMCEX) estimate of $\v$. Regarding the computational cost incurred by , inversion costs $\mathcal{O}(S^3)$, since the size of the matrix to be inverted is reduced from $NL$ to $S$. Clearly, there is no need to compute $\K_z=\K_y\otimes\K_x$. As $\S$ has binary entries, $\S\K_z\S^T$ is just a selection of $S^2$ entries in $\K_z$; and, given that $\kapz((x_i,y_j),(x_n,y_l)) =\kapx(x_i,x_n)\kapy(y_j,y_l)$, it is obtained at cost $\mathcal{O}(S^2)$. Overall, the cost incurred by is $\mathcal{O}(S^3)$. Compared to the MC approach in , the KKMCEX method is easier to implement since it only involves a matrix inversion. Moreover, since it admits a closed-form solution, it facilitates deriving bounds on the estimation error of $\hat{\v}_K$. **Remark 1**. Matrices built via the Kronecker product have been used in regression for different purposes. Related to MC, [@rao2015collaborative] leverages Kronecker product structures to efficiently solve the Sylvester equations that arise in alternating minimization iterations to find $\{\hat{\W},\hat{\H}\}$ in . On the other hand, [@pahikkala2014two; @stock2018comparative] propose a Kronecker kernel ridge regression method that can be used to extrapolate missing entries in a matrix. However, the methods in [@pahikkala2014two; @stock2018comparative] assume a complete training set and Kronecker structure for the regression matrix; this implies that the observed entries in $\bm{M}$ can be permuted to form a full submatrix. In our formulation, we introduce $\S$ which encompasses any sampling pattern in $\Omega$. Thus, the properties of the Kronecker product used in [@rao2015collaborative; @pahikkala2014two; @stock2018comparative] cannot be applied to solve $\eqref{eq:zKKMCEX}$ since $\S\K_z\S^T$ is not necessarily the Kronecker product of two smaller matrices.

**Remark 2**. The KKMCEX solution in , differs from that obtained as the solution of (\[eq:factkerw\]). On the one hand, the loss in  can be derived from the factorization-based one by using the Kronecker product kernel $\K_y\otimes\K_x$ and $\gam=\text{vec}(\B\C^T)$ to arrive at $$\begin{aligned}
& {\left|\left|P_\Omega(\M - \K_x\B\C^T\K_y)\right| \right|_\text{F}^2} \nonumber\\&  \:\: = {\left|\left|\mb - \S(\K_y\otimes\K_x)\text{vec}(\B\C^T)\right|\right|^2_2}.\end{aligned}$$ One difference between the two loss functions is that (\[eq:mcker\]) does not explicitly limit the rank of the recovered matrix $\hat{\F}=\text{unvec}({\hat{\v}_R})$ since it has $NL$ degrees of freedom through $\hat{\gam}$, while in (\[eq:factkerw\]) the rank of $\hat{\F}$ cannot exceed $p$ since $\B$ and $\C$ are of rank $p$ at most. In fact, the low-rank property is indirectly promoted in (\[eq:mcker\]) through the kernel matrices. Since $\text{rank}(\F)\leq\min(\text{rank}(\K_x),\text{rank}(\K_y))$, we can limit the rank of $\hat{\F}$ by selecting rank deficient kernels. On the other hand, the regularization terms in (\[eq:factkerw\]) and (\[eq:mcker\]) play a different role in each formulation. The regularization in (\[eq:factkerw\]) promotes smoothness on the columns of the estimated factor matrices $\{\hat{\W},\hat{\H}\}$; or, in other words, similarity between the rows of $\{\hat{\W},\hat{\H}\}$ as measured by $\kapx$ and $\kapy$. On the contrary, the regularization in (\[eq:mcker\]) promotes smoothness on $\hat{\v}$, which is tantamount to promoting similarity between the entries of $\hat{\F}$ in accordance with $\kapz$.

KKMCEX error analysis
---------------------

In order to assess the performance of KKMCEX we will rely on the mean-square error $$\label{eq:risk}
  MSE := {\mathbb{E}_{\bm e}\{ ||\v - \hat{\v}_K||^2_2\}}$$ where ${\mathbb{E}_{\bm e}\{ \cdot\}}$ denotes the expectation with respect to $\e$. Before we proceed, we will outline Nyström’s approximation.

**Definition 1**. Given a kernel matrix $\K$ and a binary sampling matrix $\S$ of appropriate dimensions, the Nyström approximation [@drineas2005nystrom] of $\K$ is $\T = \K\S^T(\S\K\S^T)^{\dagger}\S\K$, and the regularized Nyström approximation is $$\label{eq:nystreg}
  \tilde{\T} =  \K\S^T(\S\K\S^T + \mu\I)^{-1}\S\K.$$ Nyström’s approximation is employed to reduce the complexity of standard kernel regression problems such as the one in (\[eq:rralp\]). Instead of $\K$, the low-rank approximation $\T$ is used to reduce the cost of inverting large-size matrices using the MIL [@alaoui2015fast]. While it is known that the best low-rank approximation to a matrix is obtained from its top eigenvectors, Nyström’s approximation is cheaper. Using Def. 1, the following lemma provides the bias and variance of the KKMCEX estimator in (\[eq:zKKMCEX\]):

\[th:biasvar\] Given the kernel matrix $\K_z$ and its regularized Nyström approximation $\tilde{\T}_z$ with $\mu>0$, the MSE of the KKMCEX estimator is $$\begin{aligned}
\label{eq:lem1}
        \text{MSE} &= {\left|\left|(\K_z-\tilde{\T}_z)\gam\right|\right|^2_2}
        + {\mathbb{E}_{\bm e}\{ {1\over\mu^2}{\left|\left|(\K_z - \tilde{\T}_z)\S^T\bar{\e}\right|\right|^2_2}\}}
    \end{aligned}$$ where the first term accounts for the bias and the second term accounts for the variance.

Lemma \[th:biasvar\] shows that the MSE of the KKMCEX can be expressed in terms of $\tilde{\T}_z$; see proof in the Appendix. Knowing that the 2-norm satisfies ${\left|\left|\A\right|\right|^2_2}\leq{\left|\left|\A\right| \right|_\text{F}^2}$, we have $$\begin{aligned}
    &{\left|\left|(\K_z-\tilde{\T}_z)\gam\right|\right|^2_2}+ {\mathbb{E}_{\bm e}\{ {1\over\mu^2}{\left|\left|(\K_z - \tilde{\T}_z)\S^T\bar{\e}\right|\right|^2_2}\}}
    \nonumber\\&\:\:\leq {\left|\left|(\K_z-\tilde{\T}_z)\right| \right|_\text{F}^2}\left({\left|\left|\gam\right|\right|^2_2} +{\mathbb{E}_{\bm e}\{ {1\over\mu^2}{\left|\left|\S^T\bar{\e}\right|\right|^2_2}\}}\right).\end{aligned}$$ Consequently, the upper bound on the MSE is proportional to the approximation error of $\tilde{\T}_z$ to $\K_z$. This suggests selecting $\{m_{i,j}\}_{(i,j)\in\Omega}$ so that this approximation error is minimized; see also [@alaoui2015fast] where $\Omega$ is chosen according to the so-called leverage scores of $\K_z$ in order to minimize the regression error. The next theorem uses Lemma 1 to upper bound the MSE in ; see the Appendix for its proof.

\[th:mse\] Let $\sigma_{NL}$ be the maximum eigenvalue of a nonsingular $\K_z$, and $\tilde{\gam}:=\L^T\gam$, where $\L$ is the eigenvector matrix of $\K_z-\tilde{\T}_z$. If $\e$ is a zero-mean vector of iid Gaussian random variables with covariance matrix $\nu^2\I$, the MSE of the KKMCEX estimator is bounded as $$\label{eq:msebound}
MSE \leq
\frac{\mu^2\sigma_{NL}^2}{(\sigma_{NL}+\mu)^2}\sum _{i=1}^{S}\tilde{\gam_i}^2 + \sigma^2_{NL}\sum_{i=S+1}^{NL}\tilde{\gam}_i^2+\frac{S\nu^2 \sigma_{NL}^2}{\mu^2}.$$

Considering the right-hand side of (\[eq:msebound\]), the first two terms correspond to the bias, while the last term is related to the variance. We observe that when $\M$ is fully observed, that is, $S=NL$, the bias can be made arbitrarily small by having $\mu\rightarrow 0$. It is also of interest to assess how the MSE bound behaves as $S$ increases. Considering $\mu=S\mu'$ and fixed values in $(0,\infty)$ for $\mu'$, $||\tilde{\gam}_i||^2$ and $\sigma_{NL}$[^2], the bias term reduces to $$\begin{aligned}
\label{eq:biasboundtext} 
        \frac{S^2\mu'^2\sigma_{NL}^2}{(\sigma_{NL}+S\mu')^2}\sum _{i=1}^{S}\tilde{\gam_i}^2 + \sigma^2_{NL}\sum_{i=S+1}^{NL}\tilde{\gam}_i^2\;.
     \end{aligned}$$ We observe in that as $S$ increases, terms move from the second summation to the first. Therefore, whether the bias term grows or diminishes depends on the multiplication factors in front of the two summations. Since $\frac{S^2\mu'^2}{(\sigma_{NL}+S\mu')^2} \leq 1$ the bias term in decreases with $S$. On the other hand, the variance term becomes $\frac{\nu^2 \sigma_{NL}^2}{S\mu'^2}$ and decays with $S$ as well. As a result, the MSE bound in Theorem \[th:mse\] decays up until $S=NL$.

Ridge regression MCEX {#sec:RRMCEX}
=====================

Although the KKMCEX algorithm is fast when $S$ is small, the size of the matrix to be inverted in (\[eq:gamest\]) grows with $S$, hence increasing the computational cost. Available approaches to reducing the computational cost of kernel regression methods are centered around the idea of approximating the kernel matrix. For instance, [@alaoui2015fast] uses Nyström’s approximation, that our performance analysis in Section IV was based on, whereas [@yang2015randomized] relies on a sketch of $\K_z$ formed by a subset of its columns, hence reducing the number of regression coefficients; see also [@avron2017faster], where the kernel function is approximated by the inner product of random finite-dimensional feature maps, which also speeds up the matrix inversion. In this section, we reformulate the KKMCEX of Section \[sec:kmr\] to incorporate a low-rank approximation of $\K_z$ in order to obtain a reduced complexity estimate for $\bm \v$. Moreover, we also develop an online method based on this reformulation.

Recall from Eq. (\[eq:kerdot\]) that a kernel can be viewed as the inner product of vectors mapped to a feature space $\mathcal{F}_z$, namely $\kapz((x_i,y_j),(x_n,y_l)) = \langle \phi_z(x_i,y_j),\phi_z(x_n,y_l)\rangle_{\mathcal{F}_z}$. Let $\tilde{\phi}_z:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}^d$ be a feature map approximating $\kappa_z$ so that $$\label{eq:kzfeatapprox}
\kapz((x_i,y_j),(x_n,y_l)) \simeq \langle \tilde{\phi}_z(x_i,y_j),\tilde{\phi}_z(x_n,y_l)\rangle.$$ Then, we define the $NL\times d$ feature matrix $\tbphi_z:=[\tilde{\phi}_z(x_1,y_1),\allowbreak\tilde{\phi}_z(x_2,y_1),\ldots,\tilde{\phi}_z(x_N,y_L)]^T$ and form $\tK_z=\tbphi_z\tbphi^T_z$. Note that $\tK_z$ is a rank-$d$ approximation of $\K_z$, and that the equality $\K_z=\tK_z$ is only feasible when $\text{rank}(\K_z)\leq d$. Consider $\tbphi_x=[\tilde{\phi}_x(x_1),\ldots,\tilde{\phi}_x(x_N)]$ and $\tbphi_y=[\tilde{\phi}_y(y_1),\ldots,\tilde{\phi}_y(y_L)]$, where $\tilde{\phi}_x:\mathcal{X}\rightarrow \mathbb{R}^{d_x}$ and $\tilde{\phi}_y:\mathcal{Y}\rightarrow \mathbb{R}^{d_y}$, as the feature matrices forming low-rank approximations to $\K_x$ and $\K_y$, respectively. Since $\K_z=\K_y\otimes\K_x$ in KKMCEX, a prudent choice is $\tbphi_z=\tbphi_y\otimes \tbphi_x$. In the next section we will present means of constructing $\{\tbphi_x,\tbphi_y,\tbphi_z\}$ maps.

Since $\tK_z$ is a valid kernel matrix, upon replacing $\K_z$ in with $\tK_z$, the observation model reduces to

$$\bar{\m} = \S\tbphi_z\tbphi_z^T\bm\gam + \tilde{\e},$$

where $\tilde{\e} = \bar{\e} + \S(\K_z-\tK_z)\gam$. With this model, the weights in (\[eq:mcker\]) are obtained as $$\label{eq:mcfeat}
  \bm\hat{\gam} = \operatorname*{\arg\,\min}_{\gam\in \mathbb{R}^{NL}} {\left|\left|\bar{\m} - \S\tbphi_z\tbphi^T_z\gam\right|\right|^2_2} + \mu\gam^T\tbphi_z\tbphi^T_z\gam.$$ Letting $\bxi:=\tbphi_z^T\bm\gamma$ and substituting into (\[eq:mcfeat\]), we arrive at $$\label{eq:mcridge}
   \hxi=\operatorname*{\arg\,\min}_{\bxi\in \mathbb{R}^d} {\left|\left|\bar{\m} - \S\tbphi_z\bxi\right|\right|^2_2} + \mu{\left|\left|\bxi\right|\right|^2_2}$$ which admits the closed-form solution $$\label{eq:zeta}
  \hat{\bxi} = (\tbphi_z^T\S^T\S\tbphi_z+\mu\I)^{-1}\tbphi_z^T \S^T\bar{\m}.$$ Using $\hat{\bxi}$, we obtain $\hat{\v}_{R}=\tbphi_z\hat{\bxi}$ as the *ridge regression MCEX* (RRMCEX) estimate. Using the MIL  on , it follows that $$\begin{aligned}
\label{eq:zreqzk}
\hat{\bxi}=\tbphi_z^T \S^T(\S\tbphi_z\tbphi_z^T\S^T+\mu\I)^{-1}\bar{\m} \end{aligned}$$ and thus, $$\label{eq:zreqzk2}
\hat{\v}_R=\tbphi_z\hat{\bxi}=\tK_z \S^T(\S\tK_z^T\S^T+\mu\I)^{-1}\bar{\m}.$$ Therefore, shows that $\hat{\v}_R$ is equivalent to the KKMCEX solution $\hat{\v}_K$ in after replacing $\K_z$ by its low-rank approximation $\tK_z$. For error-free approximation, $\K_z=\tbphi_z\tbphi_z^T$, while $\hat{\bxi}$ in can be viewed as the primal solution to the optimization problem in , and $\hat{\gam}$ in as its dual [@shawe]. Still, obtaining $\hat{\bxi}$ requires multiplying two $d\times S$ matrices and inverting a $d\times d$ matrix, which incurs computational cost $\mathcal{O}(d^2S)$ when $S\geq d$, and $\S\tbphi_z$ is obtained at cost $\mathcal{O}(dS)$. Thus, the cost of RRMCEX grows linearly with $S$ in contrast to KKMCEX that increases with $S^3$.

By choosing an appropriate feature map so that $d\ll S$, it is possible to control the computational cost of calculating $\hat{\bxi}$. However, reduced computational cost by selecting a small $d$ might come at the price of an approximation error to $\K_z$, which correspondingly increases the estimation error of $\hat{\v}_R$. The selection of a feature matrix to minimize this error and further elaboration on the computational cost are given in Section \[sec:choosing\].

Online RRMCEX {#sec:oRRMCEX}
-------------

Online methods learn a model by processing one datum at a time. An online algorithm often results when the objective can be separated into several subfunctions, each depending on one or multiple data. In the context of MC, online implementation updates $\hat{\F}$ every time a new entry $\M_{i,j}$ becomes available. If we were to solve  each time a new observation was becoming available, inverting an $S\times S$ matrix per iteration would result in an overall prohibitively high computational cost. Still, the cost of obtaining an updated solution per observation can stay manageable using online kernel regression solvers that fall into three categories [@van2014online]: dictionary learning, recursive regression and stochastic gradient descent based. Akin to [@lu2016large; @sheikholeslami2018], we will pursue here the SGD.

Consider rewriting  entrywise as $$\label{eq:mcrr}
    \hat{\bxi}=\operatorname*{\arg\,\min}_{\bxi\in\mathbb{R}^d} \sum_{(i,j)\in \Omega}\left[m_{i,j} - \tilde{\phi}_z^T(x_i,y_j)\bxi\right]^2 + \mu{\left|\left|\bxi\right|\right|^2_2}.$$ With $n$ denoting each scalar observation, SGD iterations form a sequence of estimates $$\label{eq:sgd}
    \hat{\bxi}^n = \hat{\bxi}^{n-1} - t_n\left[\tilde{\phi}_z(x_i,y_j)(\tilde{\phi}^T_z(x_i,y_j)\hat{\bxi}^{n-1} - m_{i,j})+\mu\hat{\bxi}^{n-1}\right]$$ where $t_n$ is the step size, $n=1,\ldots,S$ and the tuple $(i,j)$ denotes the indices of the entry revealed at iteration $n$. With properly selecting $t_n$, the sequence $\hat{\bxi}^n$ will converge to  at per iteration cost $\mathcal{O}(d)$ [@bottou2012stochastic]. Apart from updating all entries in the matrix,   can also afford a simple distributed implementation using e.g., the algorithms in [@schizas].

**Remark 3**. Online algorithms for MC can be designed to solve the factorization-based formulation from (\[eq:mclag\]) rewritten as $$\begin{aligned}
\operatorname*{\arg\,\min}_{\substack{\W\in\mathbb{R}^{N\times p}\\\H\in\mathbb{R}^{N\times p}}} \!\sum_{(i,j)\in\Omega} \!\left(\!(m_{i,j} \!-\! \w_i^T\h_j)^2 \!+\! {\mu \over |\Omega^w_i|}{\left|\left|\w_i\right|\right|^2_2} \!+\! {\mu\over |\Omega^h_j|}{\left|\left|\h_j\right|\right|^2_2}\!\right)\end{aligned}$$ where $\w_i^T$ and $\h_j^T$ denote the $i^{\text{th}}$ and $j^{\text{th}}$ rows of $\H$ and $\W$ respectively, $\Omega^w_i=\{j\::\: (i,j) \in \Omega \}$, and $\Omega^h_j=\{i\::\: (i,j) \in \Omega \}$. When $m_{i,j}$ becomes available, algorithms such as SGD and online ALS update the rows $\{\w_i^T$, $\h_j^T\}$ of the coefficient matrices. This procedure can also be applied to the kernel MCEX formulation in , that solves for $\W$ and $\H$, although the rows $\{\w_i^T$, $\h_j^T\}$ cannot be updated independently due to the involvement of the kernel matrices [@zhou2012]. Then, all entries in the $i^\text{th}$ row and $j^\text{th}$ column of $\hat{\F}$ are also updated per iteration, as opposed to our method which updates the whole matrix.

Choosing the kernel matrices {#sec:choosing}
============================

In this section, we provide pointers on how to build matrices $\K_z$ for KKMCEX and $\tbphi_z$ for RRMCEX when prior information about either the matrix $\F$, or the input spaces $\mathcal{X}$ and $\mathcal{Y}$, is available.

Kernels based on the graph Laplacian {#sec:gsmc}
------------------------------------

Suppose that the columns and rows of $\F$ lie on a graph, that is, each entry of a column or row vector is associated with a node on a graph that encodes the interdependencies with entries in the same vector. Specifically, we define an undirected weighted graph $\mathcal{G}_x = (\mathcal{X}, \mathcal{E}_x, \A_x)$ for the columns of $\F$, where $\mathcal{X}$ is the set of vertices with $|\mathcal{X}|=N$, $\mathcal{E}_x \subseteq \mathcal{X}\times\mathcal{X}$ is the set of edges connecting the vertices, and $\A_x \in \mathbb{R}^{N\times N}$ is a weighted adjacency matrix. Then, functions $\{f_l: \mathcal{X} \rightarrow \mathbb{R}\}^{L}_{l=1}$ are what is recently referred to as a graph signal [@shuman2013emerging], that is, a map from the set $\mathcal{X}$ of vertices into the set of real numbers. Likewise, we define a graph $\mathcal{G}_y = (\mathcal{Y}, \mathcal{E}_y, \A_y)$ for the rows of $\F$, i.e., $\{g_n: \mathcal{Y} \rightarrow \mathbb{R}\}^{N}_{n=1}$, which are also graph signals. In a matrix of user-movie ratings for instance, we would have two graphs: one for the users and one for the movies. The graphs associated with the columns and rows yield the underlying structure of $\F$ that can be used to generate a pair of kernels.

Using $\A_x$ and $\A_y$, we can form the corresponding graph Laplacian as $\L_x:=\text{diag}(\A_x\bm 1)-\A_x$ and likewise for $\L_y$, that can serve as kernels. A family of graphical kernels results using a monotonic inverse function $r^\dagger(\cdot)$ on the Laplacian eigendecomposition as [@smola2003kernels] $$\label{eq:kerlap}
        \K = \Q r^{\dagger}(\bm\Lambda)\Q^T.$$ A possible choice of $r(\cdot)$ is the Gaussian radial basis function, which generates the diffusion kernel $r(\lambda_i)=e^{\eta\lambda_i}$, where $\lambda_i$ is the $i^\text{th}$ eigenvalue of $\L$, and $\eta$ a weight parameter. Alternatively, one can choose just the linear function $r(\lambda_i)=1 + \eta\lambda_i$, which results in the regularized Laplacian kernel. By applying different weighting functions to the eigenvalues of $\L_x$ and $\L_y$, we promote smoother or more rapidly changing functions for the columns and rows of $\hat{\F}$ [@romerospace]. While $\K_x$ and $\K_y$ are chosen as Laplacian kernels, this would not be the case for $\K_z = \K_y\otimes\K_x$ used in our KKMCEX context since it does not result from applying $r^\dagger(\cdot)$ to a Laplacian matrix. However, since $\K_x=\Q_x\Sig_x\Q_x^T$ and $\K_y=\Q_y\Sig_y\Q_y^T$, the eigendecomposition of $\K_z$ is $\K_z=(\Q_y\otimes\Q_x)(\Sig_y\otimes\Sig)(\Q^T_y \otimes\Q^T_x)$, and the notions of frequency and smoothness carry over. In other words, we are still promoting similarity among entries that are connected on the row and columns graphs through $\K_z$.

A key attribute in graph signal processing is that of “graph bandlimitedness", which arises when a signal can be generated as a linear combination of a few eigenvectors of the Laplacian matrix. Therefore, a bandlimited graph signal belongs to an RKHS that is spanned by a bandlimited kernel [@romero] that suppresses some of the frequencies of the graph. A bandlimited kernel is derived from the Laplacian matrix of a graph as in , using $$\label{eq:bandlimr}
    r(\lambda_i) = 0 \:~~ \forall i\notin\Psi,$$ where $\Psi\subseteq\mathbb{N}$ contains the indices of frequencies not to be suppressed. As mentioned earlier, we define a graph for the columns and a graph for the rows of $\F$. Therefore, graph signals contained in the columns and rows may be bandlimited with different bandwidths. In order to form $\K_z$ our KKMCEX approach, we will need to apply different weighting functions akin to the one in  for kernel matrices $\K_x$ and $\K_y$.

Kernels from known basis or features {#sec:kerfeat}
------------------------------------

In some applications the basis that spans the columns or rows of the unobserved matrix is assumed known, although this basis matrix needs not be a kernel. In order to be able to include such basis into the kernel framework, we need to generate kernel functions that span the same spaces as the columns and rows of $\F$.

Consider the input sets $\{\mathcal{X},\mathcal{Y}\}$ whose entries can be mapped into an Euclidean space through feature extraction functions $\theta_x:\mathcal{X}\rightarrow\mathbb{R}^{t_x}$ and $\theta_y:\mathcal{Y}\rightarrow\mathbb{R}^{t_y}$ such that $\theta_x(x_i):=\x_i$ and $\theta_y(y_j):=\y_j$. For instance, in a movie recommender system where the users are represented in $\mathcal{X}$ and the movies in $\mathcal{Y}$, each coordinate of $\y_j$ could denote the amount of action, drama and nudity in the movie, and $\x_i$ would contain weights according to the user’s preference for each attribute. We may then use the feature vectors to determine the similarities among entries in $\mathcal{X}$ and $\mathcal{Y}$ by means of kernel functions.

Let $\X :=[\x_1,\ldots,\x_N]^T$ and $\Y :=[\y_1,\ldots,\y_L]^T$. If $\text{span}(\F)\subseteq\text{span}(\X)$ and $\text{span}(\F^T)\subseteq\text{span}(\Y)$, we may conveniently resort to the linear kernel. The linear kernel amounts to the dot product in Euclidean spaces, which we use to define the pair $\kappa_x(x_i,x_j) = \x_i^T\x_j$ and $\kappa_y(y_i,y_j) = \y_i^T\y_j$. This leads to a straightforward construction of the kernel matrices for KKMCEX as $\K_x=\X\X^T$ and $\K_y=\Y\Y^T$.

Besides the linear kernel, it is often necessary to use a different kernel class for each $\kappa_x$ and $\kappa_y$ chosen to better fit the spaces spanned by the rows and columns of $\F$. For instance, the Gaussian kernel defined as $\kapx(x_i,x_j) = \text{exp}\{-{\left|\left|\x_i-\x_j\right|\right|^2_2}/(2\eta)\}$, is a widely used alternative in the regression of smooth functions.

Feature maps for RRMCEX
-----------------------

Aiming to construct $\tbphi_z$ that approximates $\K_z$ at reduced complexity, we choose $\tilde{\phi}_z$ with $d\ll S$. To approximate linear kernels, let $\tilde{\phi}_x(x_i)=\x_i$ and $\tilde{\phi}_y(y_j)=\y_j$ so that we can set $\tilde{\phi}_z(x_i,y_j)=\tilde{\phi}_y(y_j)\otimes\tilde{\phi}_x(x_i)$ and $\tbphi_z=\Y\otimes\X$. Note that in this case $\tbphi_z\tbphi_z^T$ yields a zero-error approximation to $\K_z=(\Y\otimes\X)(\Y\otimes\X)^T$, which renders the KKMCEX and RRMCEX solutions equivalent.

On occasion, $\X$ and $\Y$ may have large column dimension, thus rendering $\Y\otimes\X$ undesirable as a feature matrix in RRMCEX. In order to overcome this hurdle, we build an approximation to the column space of $\Y\otimes\X$ from the SVD of $\X$ and $\Y$. Consider the SVDs of matrices $\X = \U_x\D_x\V_x^T$ and $\Y = \U_y\D_y\V_y^T$, to obtain $\Y\otimes\X = (\U_y\otimes\U_x)(\D_y\otimes\D_x)(\V_y^T\otimes\V_x^T)$. Let $\tbphi_z=\U_d\D_d$, where $\U_d$ and $\D_d$ respectively hold the top $d$ singular vectors and singular values of $\Y\otimes\X$. The SVD has cost $\mathcal{O}(Nt_x^2)$ for $\X$ and $\mathcal{O}(Lt_y^2)$ for $\Y$. Comparatively, the cost of building $\K_x$ and $\K_y$ for the linear kernel is $\mathcal{O}(N^2t_x)$ and $\mathcal{O}(L^2t_y)$, respectively. Therefore, choosing RRMCEX over KKMCEX in this case incurs no extra overhead.

When a function other than the linear kernel is selected, obtaining an approximation is more complex. To approximate a Gaussian kernel on $\mathcal{X}\times\mathcal{X}$, the vectors $\{\tilde{\phi}_x(\x_i)\}^N_{i=1}$ can be obtained by means of Taylor series expansion [@cotter2011explicit] or random Fourier features [@avron2017faster], which can also approximate Laplacian, Cauchy and polynomial kernels [@avron2017faster; @rahimi2008random]. Therefore, the maps $\tilde{\phi}_x$ and $\tilde{\phi}_y$ must be designed according to the chosen kernels.

In some instances, such as when dealing with Laplacian kernels, $\X$ and $\Y$ are not available and we are only given $\K_x$ and $\K_y$. We are then unable to derive approximations to the kernel matrices by means of maps $\tilde{\phi}_x$ and $\tilde{\phi}_y$. Nevertheless, we can still derive an adequate $\tbphi_z$ to approximate $\K_z$. Indeed, Mercer’s Theorem asserts that there are eigenfunctions $\{q_n\}^{NL}_{n=1}$ in $\mathcal{H}_z$ along with a sequence of nonnegative real numbers $\{\sigma_n\}^{NL}_{n=1}$, such that $$\label{eq:mercer}
    \kapz((x_i,y_j),(x_n,y_l)) = \sum^{NL}_{n=1}\sigma_nq_n(x_i,y_j)q_n(x_n,y_l).$$ We can find from the eigendecomposition $\K_z=\bm\Q_z\Sig_z\bm\Q_z^T$, where $q_n$ is the $n\th$ eigenvector in $\Q_z$ and $\sigma_n$ the $n\th$ eigenvalue in $\Sig_z$. If $\K_z$ is low rank, we can construct $\tbphi_z=\Q_d\Sig^{{1\over 2}}_d$, where $\Q_d$ and $\Sig_d$ respectively hold the top $d$ eigenvectors and eigenvalues of $\K_z$. Note that, since $\K_z=(\Q_y\otimes\Q_x)(\Sig_y\otimes\Sig_x)(\Q^T_y \otimes\Q^T_x)$, we only need to eigendecompose smaller matrices $\K_x$ and $\K_y$ at complexity $\mathcal{O}(N^3+L^3)$. In some cases however, such as when using Laplacian kernels, the eigendecompositions are readily available, and $\bphi_z$ can be built at a markedly reduced cost.

Numerical tests
===============

In this section, we test the performance of the KKMCEX, RRMCEX and online (o)RRMCEX algorithms developed in Sections \[sec:kmr\], \[sec:RRMCEX\] and \[sec:oRRMCEX\], respectively; and further compare them to the solution of obtained with ALS [@jain2013low] and SGD [@zhou2012]. We run the tests on synthetic and real datasets, with and without noise, and measure the signal-to-noise-ratio (SNR) as $\frac{{\left|\left|\F\right| \right|_\text{F}^2}}{{\left|\left|\E\right| \right|_\text{F}^2}}$. The algorithms are run until convergence over $N_{r}=50$ realizations with different percentages of observed entries, denoted by $P_{s}=100S/(NL)$, which are taken uniformly at random per realization. As figure of merit, we use $$\text{NMSE} = \frac{1}{N_{r}}\sum\limits_{i=1}^{N_{r}}\frac{{\left|\left|\hat{\F}_i - \F\right| \right|_\text{F}^2}}{{\left|\left|\F\right| \right|_\text{F}^2}}$$ where $\hat{\F}_i$ is the estimate at realization $i$. We show results for the optimal combination of regularization and kernel parameters, found via grid search. Finally, ALS and SGD are initialized by a product of two random factor matrices, an both are stopped.

Synthetic data {#sec:synth}
--------------

We first test the algorithms on synthetic data. The $250\times250$ data matrix is generated as $\F = \K_x\bm \Gamma\K_y$, where $\bm \Gamma$ is a $250\times250$ matrix of Gaussian random deviates. For $\K_x$ and $\K_y$ we use Laplacian diffusion kernels with $\eta=1$ based on Erdös-Rényi graphs, whose binary adjacency matrices are unweighted and any two vertices are connected with probability 0.03. The resulting $\F$ is approximately low-rank, with the sum of the first 10 eigenvalues accounting for 96% of the total eigenvalue sum. Therefore, we set the rank bound $p$ to 10 for the ALS and SGD algorithms. Whether $\F$ is approximately low-rank or exactly low-rank did not affect our results, as they were similar for matrices with an exact rank of 10. For KKMCEX, $\K_z=\K_y\otimes \K_x$, and for RRMCEX $\tbphi_z=\Q_d\Sig_d^{1\over 2}$, where $\Q_d$ contains the top 250 eigenvectors of $\K_z$, and $\Sig_d$ the corresponding top 250 eigenvalues.

Fig. \[fig:erd\] shows the simulated NMSE when $\M$ is noiseless (a) or noisy (b). We deduce from Fig. \[fig:erd\]a that all algorithms except SGD achieve a very small NMSE, below 0.003 at $P_s=1\%$ that falls to 0.007 at $P_s=10\%$. Of the three algorithms, KKMCEX has the smallest error except at $P_s=1\%$, where RRMCEX performs best. Although the error drops below 0.005 for SGD at $P_s > 4\%$, it is outperformed by the other algorithms by an order of magnitude. Fig. \[fig:erd\]b shows the same results when Gaussian noise is added to $\F$ at $\text{snr}=1$. We observe that KKMCEX and RRMCEX are matched and attain the lowest error, whereas ALS and SGD have larger errors across $P_s$. This corroborates that thanks to the regularization term that smoothes over all the entries instead of row or column-wise, the noise effect is reduced. Interestingly, RRMCEX is able to reduce the noise effect despite the bias it suffers because it only uses the top 250 eigenvalues of $\K_z$ from a total of $62{,}500$. This is mainly due to the additive noise being evenly distributed across the eigenspace of $\K_z$. Therefore, by keeping only the eigenvectors associated with the top 250 eigenvalues in $\K_z$, we are discarding those dimensions in which the SNR is lower.

Fig. \[fig:erdtime\] depicts the time needed for the algorithms to perform the simulations reported in Fig. \[fig:erd\]. We observe in Fig. \[fig:erdtime\]a that RRMCEX has an almost constant computation time, whereas the time for KKMCEX grows with $P_s$ as expected since the size of the matrix to be inverted increases with $S$. On the other hand, ALS and SGD require less time than KKMCEX for the larger values of $P_s$, but are always outperformed by RRMCEX. Moreover, the ALS time is reduced as $P_s$ increases because the number of iterations required to converge to the minimum is smaller. Fig. \[fig:erdtime\]b suggests that the noise only impacts ALS, which has its computation time rise considerably across all $P_s$. Overall, Figures \[fig:erd\] and \[fig:erdtime\] illustrate that RRMCEX has the best performance for the synthetic matrix both in terms of NMSE and computational cost.

Temperature measurements {#sec:real}
------------------------

In this case, $\F$ has size $150\times365$ comprising temperature readings taken by 150 stations over 365 days in 2002 in the United States[^3]. The columns and rows of $\F$ are modeled as graph signals with $\A_x$ and $\A_y$ for the graphs formed by the stations and the days of the year, respectively. We use Laplacian diffusion kernels both for $\K_x$ and $\K_y$, while $\K_z$ and $\tbphi_z$ are obtained as in the tests on synthetic data, except that $\tbphi_z$ is constructed with the top 150 eigenvectors of $\K_z$. The matrix $\A_x$ is obtained as in [@chen], where a graph $\mathcal{G}$ with unweighted adjacency matrix $\P$ is generated for the stations, and each station is a vertex connected to the 8 geographically closest stations. Next, we obtain the undirected graph $\mathcal{G}'$ with symmetric adjacency matrix $\P’ =\text{sign}(\P^T+\P)$. Finally, the entries of $\A_x$ are constructed as $(\A_x)_{i,j}= \text{exp}(-\frac{N^2 d_{i,j}}{\sum_{i,j} d_{i,j}})$, where $\{d_{i,j}\}$ are geodesic distances on $\mathcal{G}$. We adopt a graph on which each day is a vertex and each day is connected to the 10 past and future days to form $\A_y$.

Fig. \[fig:time\] shows the simulated tests for (a) the matrix of temperature readings, and (b) the same matrix with additive Gaussian noise at $\text{snr}=1$. Fig. \[fig:time\]a demonstrates that KKMCEX achieves the lowest error for the first three data points, while afterwards ALS has a slight edge over KKMCEX. The real data matrix $\F$ is approximately low rank, since the sum of the first 10 singular values accounts for 75% of the total sum. This explains why RRMCEX fares worse than KKMCEX. Because it only contains the top 150 eigenvectors of $\K_z$, which is full rank, the vectorized data $\m$ lies in part outside the space spanned by $\tbphi_z$. Indeed, increasing the number of eigenvectors in $\tbphi_z$ results in a lower error, although the computational cost increases accordingly. Fig. \[fig:time\]b further demonstrates that the addition of noise has the least impact on RRMCEX, which attains the lowest error slightly below KKMCEX. On the other hand, ALS has a marginally higher error whereas the gap between SGD and the other three methods remains. Fig. \[fig:timetime\] depicts the computational time for the results in Fig. \[fig:time\]a, which are similar to those obtained for the synthetic dataset.

Mushroom dataset
----------------

The Mushroom dataset[^4] comprises 8,124 labels and as many feature vectors. Each label indicates whether a sample is edible or poisonous, and each vector has 22 entries describing the shape, color, etc. of the mushroom sample. After removing items with missing features, we are left with 5,643 labels and feature vectors. Here, we solve a clustering problem in which $\F$ is a $5,643\times 5,643$ adjacency matrix, where $\F_{i,j} = 1$ if the $i^{\text{th}}$ and the $j^{\text{th}}$ mushroom samples belong to the same class (poisonous or edible), and $\F_{i,j} = -1$ otherwise.

We encode the matrix stacking the feature vectors via one-hot encoding to produce a $5,643\times 98$ binary feature matrix analogous to $\X$ in Section \[sec:kerfeat\]. We build the kernel matrix $\K_x$ from the Pearson correlation coefficients of the rows of $\X$, and let $\K_y=\K_x$. The feature matrix $\tbphi_z$ for RRMCEX is built using the top 3,000 left singular vectors of $\X\otimes\X$.

Fig. \[fig:mushroom\]a shows the test results on the mushroom adjacency matrix from $S=2,000$ $(P_s=0.006\%)$ to $S=20{,}000$ $(P_s=0.036\%)$ in steps of 1,000 observations. KKMCEX and RRMCEX achieve similar NMSE, while SGD has an error one order of magnitude higher, and ALS outperforms both by around one order of magnitude. This difference with ALS is because regression-based methods restrict the solution to belong to the space spanned by the basis matrix. On the other hand, when solving , we do not enforce the constraints $\W\in\mathcal{H}_x$, $\H\in\mathcal{H}_y$ [@bazerque2013; @zhou2012]. Therefore, when the prior information encoded in the kernel matrices is imperfect, ALS might be able to find a factorization that fits better the data at the cost of having $\hat{\W}\notin\mathcal{H}_x$ and $\hat{\H}\notin\mathcal{H}_y$. However, in Fig. \[fig:mushroom\]b we see that the computational cost for ALS and SGD is much higher than for KKMCEX and RRMCEX for the smaller $\P_s$. On the other hand, the time for ALS decreases with $S$ due to requiring les iterations to converge, whereas for KKRRMCEX and RRMCEX it increases with $S$.

Online MC
---------

In the online scenario, we compare the (o)RRMCEX algorithm with online (o)ALS and SGD. One observation is revealed per iteration at random, and all three algorithms process a single observation per iteration in a circular fashion. Per realization, we run tests on both synthetic and temperature matrices with $P_s=10\%$, that is, $S=6{,}250$ and $S=5{,}475$ observations for the synthetic and temperature matrices, respectively, for a single realization.

Fig. \[fig:sg\]a depicts the tests for the noiseless synthetic matrix. Clearly, (o)RRMCEX converges much faster than SGD and (o)ALS. Indeed, as opposed to SGD and (o)ALS, which require several passes over the data, (o)RRMCEX approaches the minimum in around $6{,}000$ iterations. Moreover, it achieves the smallest NMSE of 0.0004, which is slightly below the 0.0011 obtained by SGD. Fig. \[fig:sg\]b shows the results for the temperature matrix without noise. Again, we observe that (o)RRMCEX converges the fastest to the minimum, whereas SGD requires many passes through the data before it starts descending, while (o)ALS converges much faster than with the synthetic data. Regarding the NMSE, (o)RRMCEX and SGD achieve the same minimum value.

The tests on the Mushroom dataset are run with $S=10{,}000$ $(P_s=0.033\%)$ and $S=20{,}000$ $(P_s=0.036\%)$ observations following the same procedure as with the synthetic and temperature datasets. Fig. \[fig:mushonline\] shows results for the Mushroom adjacency matrix with the error for $S=20{,}000$ plotted in solid lines, and for $S=10{,}000$ in dotted lines. We observe that for $S=20{,}000$, (o)RRMCEX crosses the minimum of (o)ALS and SGD in 7 seconds, whereas (o)ALS and SGD converge to this minimum in 12 and 200 seconds, respectively. Afterwards, the line for (o)RRMCEX keeps descending until an error of $0.012$ is reached. For $S=10{,}000$ the convergence time of (o)RRMCEX and SGD remains almost unchanged, whereas for (o)ALS it increases to 26 seconds. Moreover, the error of both (o)ALS and SGD grows much larger, whereas (o)RRMCEX exhibits just a small increase.

Conclusions
===========

In this paper, we have taken a comprehensive look at MC under the framework of RKHS. We have viewed columns and rows of the data matrix as functions from an RKHS, and leveraged kernel theory to account for the available prior information on the contents of the sought matrix. Moreover, we have developed two estimation algorithms that offer simplicity and speed as their main advantages. When the number of observed data is small, KKMCEX obtains the full matrix estimate by inverting a reduced-size matrix thanks to the Representer Theorem. On the other hand, when the number of observations is too large for KKMCEX to handle, RRMCEX can be employed instead in order to lower the computational cost with no impact on the recovery error when noise is present. In addition, RRMCEX can be easily turned into an online method implemented via SGD iterations. Compared to mainstream methods designed for the factorization-based formulation, namely ALS and SGD, our RRMCEX exhibited improved performance in simulated and real data sets.

Our future research agenda includes improving both KKMCEX and RRMCEX through parallel and accelerated regression methods, as well as designing robust sampling strategies for MCEX formulated as a kernel regression.

Proof of Lemma 1
----------------

For the KKMCEX estimator (\[eq:zKKMCEX\]), the MSE is given as $$\label{eq:risk}
  MSE := {\mathbb{E}_{\bm e}\{ {\left|\left|\v - \hz_K\right|\right|^2_2}\}} = {\mathbb{E}_{\bm e}\{ {\left|\left|\v - \K_z\hgam\right|\right|^2_2}\}}.$$ Plugging the estimator from (\[eq:gamest\]) into (\[eq:risk\]) yields $$\begin{aligned}
  MSE &= {\mathbb{E}_{\bm e}\{ {\left|\left|\v-\K_z \S^T (\S\K_z\S^T + \mu\I)^{-1}(\S\v+\bar{\e})\right|\right|^2_2}\}}\nonumber \\
  &= {\left|\left|(\I - \K_z \S^T (\S\K_z\S^T + \mu\I )^{-1}\S)\v\right|\right|^2_2}\nonumber \\ &+ {\mathbb{E}_{\bm e}\{ {\left|\left|\K_z\S^T(\S\K_z\S^T+\mu\I)^{-1}\bar{\e}\right|\right|^2_2}\}} \label{eq:biasvar}\end{aligned}$$ where we have used that $\mathbb{E}\{\e\}=\bm 0$. Further, the first and second terms in (\[eq:biasvar\]) are the bias and variance of the KKMCEX estimator, respectively. If we substitute $\v=\K_z\gam$ into the first term of (\[eq:biasvar\]), we obtain $$\begin{aligned}
bias &= {\left|\left|(\I - \K_z\S^T (\S\K_z\S^T + \mu\I)^{-1}\S )\K_z\gam\right|\right|^2_2}\nonumber \\
&= {\left|\left|(\K_z - \K_z \S^T (\S\K_z\S^T + \mu\I)^{-1}\S\K_z )\gam\right|\right|^2_2}\nonumber \\
& = {\left|\left|(\K_z-\tilde{\T}_z)\gam\right|\right|^2_2}\label{eq:bias}\end{aligned}$$ where $\tilde{\T}_z$ is the regularized Nyström approximation of $\K_z$ in (\[eq:nystreg\]). On the other hand, the variance term is $$\begin{aligned}
var &= {\mathbb{E}_{\bm e}\{ {\left|\left|\K_z\S^T(\S\K_z\S^T+\mu\I)^{-1}\bar{\e}\right|\right|^2_2}\}}\nonumber \\
&\mspace{-25mu}= \mathbb{E}\{{1\over\mu^2}\left|\left|\K_z\S^T(\S\K_z\S^T+\mu\I)^{-1}\right.\right.\nonumber\\&\left.\left. \:\:\:(\mu\I+\S\K_z\S^T-\S\K_z\S^T)\bar{\e}\right|\right|_2^2\}\nonumber \\
&\mspace{-25mu}= {\mathbb{E}_{\bm e}\{ {1\over\mu^2}{\left|\left|\K_z\S^T - \K_z\S^T(\S\K_z\S^T+\mu\I)^{-1}\S\K_z\S^T\bar{\e}\right|\right|^2_2}\}}\nonumber\\
&\mspace{-25mu}= {\mathbb{E}_{\bm e}\{ {1\over\mu^2}{\left|\left|(\K_z - \K_z\S^T(\S\K_z\S^T+\mu\I)^{-1}\S\K_z)\S^T\bar{\e}\right|\right|^2_2}\}}\nonumber\\
&\mspace{-25mu}= {\mathbb{E}_{\bm e}\{ {1\over\mu^2}{\left|\left|(\K_z - \tilde{\T}_z)\S^T\bar{\e}\right|\right|^2_2}\}}.\label{eq:quad}\end{aligned}$$ Adding the two terms in  and , we obtain the MSE in .

Proof of Theorem 2
------------------

Since $\K_z - \tilde{\T}_z$ appears in the bias and variance terms in Lemma 1, we will first derive an upper bound on its eigenvalues that will eventually lead us to a bound on the MSE. To this end, we will need a couple of lemmas.

\[lem:sim\] Given a symmetric matrix $\A\in\mathbb{R}^{N\times N}$ and a symmetric nonsingular matrix $\B\in\mathbb{R}^{N\times N}$, it holds that $\lambda_k(\A\B) = \lambda_k(\B^{1\over 2}\A\B^{1\over 2})$; and also $\lambda_k(\A\B) \leq \lambda_k(\A)\lambda_N(\B)$.

Since $\B$ is invertible and symmetric, we can write $\A\B = \B^{-{1\over 2}}(\B^{1\over 2}\A\B^{1\over 2})\B^{1\over 2}$. Therefore, $\A\B$ is similar to $\B^{{1\over 2}}\A\B^{1\over 2}$, and they both share the same eigenvalues. Let $\mathcal{U}\subset\mathbb{R}^{N}\setminus\{\bm 0\}$. From the min-max theorem [@lax], the $k^{th}$ eigenvalue of $\A$ satisfies $$\lambda_k(\A) = \min_{\mathcal{U}} \left\{ \max_{\x\in \mathcal{U}} \frac{\x^T\A\x}{\x^T\x} \:|\: \text{dim}(\mathcal{U})=k \right\}.$$ Therefore, we have $$\begin{aligned}
    \lambda_k(\A\B) &= \lambda_k(\B^{1\over 2}\A\B^{1\over 2}) \nonumber\\&\mspace{-50mu}=
      \min_{\mathcal{U}} \left\{ \max_{\x\in \mathcal{U}} \frac{\x^T\B^{1\over 2}\A\B^{1\over 2}\x}{\x^T\x} \:|\: \text{dim}(\mathcal{U})=k \right\} \nonumber\\
    &\mspace{-50mu}= \min_{\mathcal{U}} \left\{ \max_{\x\in \mathcal{U}} \frac{\x^T\B^{1\over 2}\A\B^{1\over 2}\x}{\x^T\B^{1\over 2}\B^{1\over 2}\x} \frac{\x^T\B\x}{\x^T\x} \:|\: \text{dim}(\mathcal{U})=k \right\}\nonumber \\
    &\mspace{-50mu}\leq \min_{\mathcal{U}} \left\{ \max_{\x\in \mathcal{U}} \frac{\x^T\A\x}{\x^T\x} \:|\: \text{dim}(\mathcal{U})=k \right\}\lambda_N(\B)\nonumber \\
    &\mspace{-50mu}= \lambda_k(\A)\lambda_N(\B).\end{aligned}$$

The following lemma bounds the eigenvalues of $\K_z - \tilde{\T}_z$, and the regularized Nystrom approximation $\tilde{\T}_z$ in (\[eq:nystreg\]).

\[lem:nyst\] With $\K_z$ as in  and $\tilde{\T}_z$ as in , the eigenvalues of $\K_z-\tilde{\T}_z$ are bounded as $$\K_z - \tilde{\T}_z \preceq \frac{\mu\sigma_{NL}}{\sigma_{NL}+\mu}\I_S' + \sigma_{NL}\I_S$$ where $\sigma_{NL}$ is the largest eigenvalue of $\K_z$, $\I_S := \text{diag}([0,\allowbreak 0,\ldots,1,1])$ has $S$ zeros on its diagonal, and $\I'_S := \I-\I_S$.

Using the eigendecomposition $\K_z=\Q_z\Sig_z\Q_z^T$, we can write $$\begin{aligned}
    \K_z-\tilde{\T}_z &= \K_z - \K_z\S^T(\S\K_z\S^T+\mu\I)^{-1}\S\K_z \nonumber \\
    &= \Q_z\Sig_z^{1\over 2}\left[\I-\Sig_z^{1\over 2}
    \Q_z^T\S^T(\S\Q_z\Sig_z^{1\over 2}\Sig_z^{1\over 2}
    \Q_z^T\S^T \right.\nonumber \\&\left.\:\:\:\:+\mu \I)^{-1}\S\Q_z\Sig_z^{1\over 2}\right]
    \Sig_z^{1\over 2}\Q_z^T.\label{eq:midmat}\end{aligned}$$ Applying the MIL to the matrix inside the square brackets of , we arrive at $$\begin{aligned}
    &\I-\Sig_z^{1\over 2}
    \Q_z^T\S^T(\S\Q_z\Sig_z^{1\over 2}\Sig_z^{1\over 2}
    \Q_z^T\S^T+\mu \I)^{-1}\S\Q_z\Sig_z^{1\over 2} \nonumber\\&\:\:= (\I + {1\over\mu}\Sig_z^{1\over 2}\Q_z^T\S^T\S\Q_z\Sig_z^{1\over 2})^{-1}.\end{aligned}$$ That in turn implies $$\begin{aligned}
\label{eq:kzp}
    &\K_z-\tilde{\T}_z = \mu\Q_z\Sig_z^{1\over 2}(\mu\I + \Sig_z^{1\over 2}\Q_z^T\S^T\S\Q_z\Sig_z^{1\over 2})^{-1}\Sig_z^{1\over 2}\Q_z^T \nonumber \\
    &= \mu\Q_z\Sig_z^{1\over 2}(\Sig_z + \mu\I  - \Sig_z + \Sig_z^{1\over 2}\Q_z^T\S^T\S\Q_z\Sig_z^{1\over 2})^{-1}\Sig_z^{1\over 2}\Q_z^T \nonumber \\
  &= \mu\Q_z\Sig_z^{1\over 2}\left[(\Sig_z+\mu\I)^{1\over2}\left(\I -  (\Sig_z+\mu\I)^{-{1\over2}}\Sig_z(\Sig_z+\mu\I)^{-{1\over 2}}\right.\right.\nonumber\\& \left.\left. \:\:\:\:+ (\Sig_z+\mu\I)^{-{ 1\over 2}}\Sig_z^{1\over 2}\Q_z^T\S^T\S\Q_z\Sig_z^{1\over 2}(\Sig_z+\mu\I)^{-{1\over 2}})\right)\right.\nonumber \\&\left. \:\:\:\:(\Sig_z+\mu\I)^{1\over2}\right]^{-1}\Sig_z^{1\over 2}\Q_z^T\nonumber \\
  &= \mu\Q_z\Sig_z^{1\over 2}(\Sig_z+ \mu\I)^{-{1\over2}}(\I-\P)^{-1}(\Sig_z+ \mu\I)^{- {1\over2}}\Sig_z^{1\over 2}\Q_z^T \end{aligned}$$ where $$\begin{aligned}
\label{eq:P}
    \P =&\Sig_z(\Sig_z+\mu\I)^{-1} \\&- (\Sig_z+\mu\I)^{-{ 1\over 2}}\Sig^{1\over 2}\Q_z^T\S^T\S\Q_z\Sig^{1\over 2}(\Sig_z+\mu\I)^{-{1\over 2}}.\nonumber \end{aligned}$$ Regarding the eigenvalues of $\K_z-\tilde{\T}_z$ in , we can bound them as $$\begin{aligned}
    &\lambda(\K_z-\tilde{\T}_z) \nonumber\\&= \mu\,\lambda(\Q_z\Sig_z^{1\over 2}(\Sig_z+ \mu\I)^{-{1\over2}}(\I-\P)^{-1}(\Sig_z+ \mu\I)^{- {1\over2}}\Sig_z^{1\over 2}\Q_z^T)\nonumber\\
     &= \mu\,\lambda(\Sig_z^{1\over 2}(\Sig_z+ \mu\I)^{-{1\over2}}(\I-\P)^{-1}(\Sig_z+ \mu\I)^{- {1\over2}}\Sig_z^{1\over 2})\nonumber\\
     &= \mu\,\lambda((\I-\P)^{-1}(\Sig_z+ \mu\I)^{-1}\Sig_z)\nonumber\\
     &\leq\frac{\mu\sigma_{NL}}{\sigma_{NL}+\mu}\lambda((\I-\P)^{-1})\label{eq:kzbound1}
 \end{aligned}$$ where $\lambda(\cdot)$ denotes the eigenvalues of a matrix, and we have applied Lemma \[lem:sim\] on the third equality and the last inequality. Knowing that $\lambda(\I-\P)=\I-\lambda(\P)$ we can now bound the eigenvalues of $\P$ as $$\begin{aligned}
\label{eq:C2}
    \lambda(\P) &\!=\!\lambda(\Sig_z^{1\over 2}(\Sig_z\!+\!\mu\I)^{- {1\over 2}}\Q_z^T(\I \!-\! \S^T\S)\Q_z(\Sig_z\!+\!\mu\I)^{- {1\over 2}}\Sig_z^{1\over 2})\nonumber\\
 &=\lambda(\Q_z^T(\I - \S^T\S)\Q_z(\Sig_z+\mu\I)^{- {1}}\Sig_z)\nonumber\\
 &\leq \frac{\sigma_{NL}}{\sigma_{NL}+\mu} \lambda(\Q_z^T(\I - \S^T\S)\Q_z)\nonumber\\
 &= \frac{\sigma_{NL}}{\sigma_{NL}+\mu}\lambda(\I_S)\end{aligned}$$ where we have applied Lemma \[lem:sim\] on the second and third inequalities, and $\I_S:=\text{diag}{[0,0,\ldots,1,1]}$ has $S$ zeros on its diagonal. Next, we have that $\I-\P\succeq \I - \frac{\sigma_{NL}}{\sigma_{NL}+\mu}\I_S$, and thus $$\label{eq:iqbound}
(\I-\P)^{-1}\preceq  \I_S' + \frac{\sigma_{NL}+\mu}{\mu}\I_S$$ where $\I_S':=\I-\I_S$. Finally, combining  with  yields $$\label{eq:kt}
    \K_z - \tilde{\T}_z \preceq \frac{\mu\sigma_{NL}}{\sigma_{NL}+\mu}\I_S' + \sigma_{NL}\I_S$$ which concludes the proof.

Using Lemmas \[lem:sim\] and \[lem:nyst\], we can proceed to establish a bound on the bias and variance. Considering the eigendecomposition $\K_z-\tilde{\T}_z=\L\bm\Lambda\L^T$, we can write the bias in (\[eq:bias\]) as $$\begin{aligned}
bias &= {\left|\left|\L\bm\Lambda\L^T\gam\right|\right|^2_2}.\end{aligned}$$ With $\tilde{\gam}:=\L^T\gam$, and using Lemma \[lem:nyst\] the bias is bounded as $$\begin{aligned}
bias
  &  = {\left|\left|\L\bm\Lambda\tilde{\gam}\right|\right|^2_2}
   = \tilde{\gam}^T\bm\Lambda^2\tilde{\gam}\nonumber\\
&\leq \frac{\mu^2\sigma_{NL}^2}{(\sigma_{NL}+\mu)^2}\tilde{\gam}^T\I_S'\tilde{\gam} + \sigma_{NL}^2\tilde{\gam}^T\I_S\tilde{\gam}\nonumber\\
&= \frac{\mu^2\sigma_{NL}^2}{(\sigma_{NL}+\mu)^2}\sum _{i=1}^{S}\tilde{\gam_i}^2 + \sigma^2_{NL}\sum_{i=S+1}^{NL}\tilde{\gam}_i^2.\label{eq:biasbound}\end{aligned}$$

To bound the variance in (\[eq:quad\]), recall that $\e$ is a Gaussian random vector with covariance matrix $\nu^2\I$, while $\bar{\e}$ has covariance matrix $\nu^2\S\S^T$. Then (\[eq:quad\]) is a quadratic form in $\e$, whose the variance becomes $$\begin{aligned}
\label{eq:var1}
    var &={\mathbb{E}_{\bm e}\{ {1\over\mu^2}{\left|\left|(\K_z - \tilde{\T}_z)\S^T\bar{\e}\right|\right|^2_2}\}} \nonumber\\&= {\nu^2\over \mu^2}{\text{Tr}(\S(\K_z - \tilde{\T}_z)^2\S^T)}\nonumber\\&= {\nu^2\over \mu^2}{\text{Tr}((\K_z - \tilde{\T}_z)^2\S^T\S)}.\end{aligned}$$ The matrix inside the trace in  has $NL-S$ zero entries in its diagonal. Lemma \[lem:nyst\], on the other hand, implies that diagonal entries of $\K_z-\tilde{\T}_z$ are smaller than its largest eigenvalue; that is, $\left[\K_z-\tilde{\T}_z\right]_{i,i}\leq\sigma_{NL}$. Coupling this with  yields $$\begin{aligned}
var
        \leq \frac{S\nu^2 \sigma_{NL}^2}{\mu^2}.\label{eq:varbound}\end{aligned}$$ Finally, combining the bias bound in (\[eq:biasbound\]) with the variance bound in (\[eq:varbound\]), yields the bound for the MSE as $$MSE \leq
\frac{\mu^2\sigma_{NL}^2}{(\sigma_{NL}+\mu)^2}\sum _{i=1}^{S}\tilde{\gam_i}^2 + \sigma^2_{NL}\sum_{i=S+1}^{NL}\tilde{\gam}_i^2+\frac{S\nu^2 \sigma_{NL}^2}{\mu^2}.$$

[^1]: This work was supported by the Ministerio de Economia y Competitividad of the Spanish Government and ERDF funds (TEC2016-75067-C4-2-R,TEC2015-515 69648-REDC), Catalan Government funds (2017 SGR 578 AGAUR), and NSF grants (1500713, 1514056, 1711471 and 1509040).

[^2]: Note that $||\tilde{\gam}_i||^2$ and $\sigma_{NL}$ depend on the selected kernel $\K_z$ and matrix $\F$, and do not depend on $\S$.

[^3]: http://earthpy.org/ulmo.html

[^4]: http://archive.ics.uci.edu/ml
