text:
\section{Introduction}
Dimensionality reduction is often an essential step before applying a machine learning method. It can be equivalently defined as an embedding of the data in a lower dimensional space by preserving certain important features.
Diffusion maps~\cite{DiffusionMaps,Coifman05geometricdiffusions} provides a meaningful embedding of high dimensional data since they are constructed to preserve an isometry relation. A central object involved in the definition of this isometry is a diffusion kernel defining a diffusion process on the points of the dataset.
Given $N$ data points, a diffusion embedding in $\mathbb{R}^N$ is obtained thanks to the spectral decomposition of the diffusion kernel. Then, the eigenvalue decomposition is truncated in order to obtain an approximate embedding in a lower dimensional Euclidean space. For an $m$-dimensional embedding space, spectral truncation is the optimal rank-$m$ approximation of the diffusion kernel in certain sense as it is explained in~\cite{Nadler}. In the same context, given a cloud of data points in a measurable space, we propose to obtain another positive semi-definite kernel matrix by solving a semi-definite program (SDP) involving the same diffusion kernel generating the diffusion process.
The solution of this optimization problem is often of low rank as motivated hereafter.

Coming back to the embedding problem described before, we should emphasize that the SDP embedding method proposed in this paper naturally shares similarities with Diffusion Maps and can be categorized as a non-linear dimensionality reduction method yielding a data dependent kernel as a natural consequence. Let us mention some related works in the literature.
The motivation of our work can be related to recent papers concerning SDP and estimation problems on graphs.
Indeed, several authors studied the connection between an SDP and the max-cut problem, $\mathbb{Z}_2$ or angular synchronization, as well as community detection \cite{Burer2003,BoumalGeNPowMeth,MontanariPNAS}. Also, in the context of the angular synchronization problem, the paper~\cite{Bandeira2016} investigates when a SDP relaxation has a global rank one solution which solves exactly an angular synchronization problem. Incidentally, we may also mention that an approximation of a $k$-means clustering can also be formulated as an SDP~\cite{PengWei}. The conclusions presented here can be viewed as a generalization of previous ideas applied to graph data~\cite{Burer2003,BoumalGeNPowMeth,Bandeira2016,MontanariPNAS} in the context of kernel methods.

Before going to the formal exposition of the proposed method, let us first summarize our proposal from a ``practitioner'' viewpoint. Three major aspects of our methods can be outlined as follows:
\begin{itemize}
\item From a practical perspective, we are given a set of $N$ data points $\Omega=\{x_1,\dots, x_N\}\subset\mathbb{R}^d$ that are going to be embedded in a low dimensional vector space, say in $\mathbb{R}^r$. Typically, we assume $r\leq d$. The dimension $r<N$ of this space is determined by the numerical solution of the convex problem defined hereafter. Associated with this data, a positive semi-definite (\emph{p.s.d.}) kernel matrix $K\in \mathbb{R}^{N\times N}$ has to be chosen, an interesting choice being a kernel related to a diffusion process (see section~\ref{sec:Diff}). Then, we phrase a semi-definite program for finding a unique optimal \emph{p.s.d.} kernel matrix $\rho^\star\in \mathbb{R}^{N\times N}$.  The ``SDP embedding'' of the data $\Xi: \mathbb{R}^d \to \mathbb{R}^r$ is then given by the eigendecomposition $\rho^\star = \Xi\Xi^\intercal$, where the dimension $r$ of the embedding is merely the rank of $\rho^\star$.
 Let us comment about the embedding provided by $\rho^\star$ and the kernel distance associated to it. Indeed, the squared distance between two $x,y\in \Omega$ in the embedding space is given by
\[
d^2_{\rho^\star}(x,y) \triangleq \|\Xi(x)-\Xi(y)\|_2^2 = \rho^\star(x,x) + \rho^\star(y,y) -2\rho^\star(x,y).
\]
In order to ``learn'' $\rho^\star$, we propose to solve the following optimization problem
\[
\min_{\rho\succeq 0}\sum_{x,y\in \Omega}K(x,y)d^2_{\rho}(x,y)\ {\rm subject\ to}\ \rho(x,x) =K(x,x)\ {\rm for\ all}\ x\in \Omega.
\]
This ``training'' problem yields a solution $\rho^\star$ such that the distance $d_{\rho^\star}(x,y)$  is small in the embedding space if $K(x,y)>0$ and is large if $K(x,y)<0$, i.e, if $x$ and $y$ are not similar. Noticeably, the embedding dimension is given directly by the rank of $\rho^\star$ and it is not determined in advance as it is often the case for dimensionality reduction methods.
As explained hereafter, this semi-definite program may be equivalently formulated as a nuclear norm minimization problem subject to linear constraints which motivates why empirically the optimal solution $\rho^\star\in \mathbb{R}^{N\times N}$ is often a very low rank matrix. More specifically, the optimization problem promotes solutions of low rank and, therefore, a low dimensional embedding.
\item From the computational viewpoint, the semi-definite program is not solved directly, although it is possible for moderate scale problems using interior point methods~\cite{Nesterov}.  Specifically, the numerical method used in this paper solves a rank constrained version based on the factorization of $\rho^\star = HH^\intercal$ with $H\in \mathbb{R}^{N\times r_0}$.  The calculation of the embedding coordinates includes two steps. An  initial seed is chosen as a random $N\times r_0$ matrix with $r_0$ taken in practice such that $5\leq r_0< 20$. Then, a projected power method as proposed in~\cite{BoumalGeNPowMeth} is applied until convergence. In practice, the resulting matrix has often a rank $r \leq 5$. A dual certificate is used in order to verify the optimality of the candidate solution. Secondly, a singular value decomposition of the $N\times r_0$ matrix obtained is performed in order to obtain the embedding coordinates $\Xi$.
\item In practice, after this initial ``training'' phase on $N$ data points, in the ``extension'' phase, additional points can be embedded with a low computational cost. Indeed, the embedding coordinates $\Xi$ admit a natural out-of-sample extension and yield an extension of the kernel matrix $\rho$ to a kernel function by relying on the conditions for optimality given in terms of a dual certificate. Concretely, the kernel matrix $\rho^\star\in\mathbb{R}^{N\times N}$ can be extended to a kernel function $\bar{\rho}^\star:\mathbb{R}^d\times \mathbb{R}^d\to \mathbb{R}$. In Proposition~\ref{Prop:Extension}, we argue that the extension of the kernel matrix is very natural. Concerning the extension formula for the embedding coordinates, it differs from the previous approaches based on Nystr\"om-type methods~\cite{WilliamSeeger} by its non-linear character. Hence, we call the extension a ``projected Nystr\"om extension''.
\end{itemize}
The properties of the embedding method are twofold.
Firstly, the length of the embedding vectors is fixed as a constraint, i.e., the data points are always embedded within a  spherical shell determined in advance. However, the dimension of the embedding space remains unconstrained in theory. Furthermore, due to the choice of the diffusion kernel, the embedding has a stratigraphic structure, i.e., data points corresponding to the same probability at stationarity (w.r.t. the diffusion process) are positioned on the same sphere (for an illustration see~Figure~\ref{Fig:OutliersBandwidth}). Hence, we observe empirically that our method is relatively robust to the presence of outliers, as it can be expected from other works involving similar SDP's~\cite{MontanariPNAS}. Indeed, in contrast to  the SDP approach, it was observed in the context of community detection in sparse graphs that spectral methods often suffer from the localization of the eigenvectors~\cite{MontanariPNAS}.
Secondly, the embedding coordinates satisfy a ``mean value property'' (see section~\ref{sec:Diff}) motivating the convenient out-of-sample formula for any new points.
In a more general context, we also discuss in this paper a more general optimization problem: finding a Mercer kernel function $k(x,y)$, rather than a Mercer kernel matrix (see, for instance,~\cite{LanckrietLearningKernel} for another approach to kernel  matrix learning.). Let us briefly motivate this question. Indeed, kernel based methods require the choice a kernel function $k(x,y)$, an example being the framework of function estimation in a reproducing kernel Hilbert space. Then, the training problem is defined by using a kernel (Gram) matrix $[k(x_i,x_j)]_{i,j}$, while the solution of the estimation procedure is often of the form $f(\cdot) = \sum_{i =1}^n \alpha_i k(x_i,\cdot)$, with $\alpha_i\in \mathbb{R}$. Hence, the selection of an optimal kernel function is relevant in this context.

Let us briefly outline the structure of this paper in order to provide a reading guide.
\begin{itemize}
\item[-] In Section~\ref{sec:2}, an infinite dimensional analogue of the SDP studied in this paper is discussed, in order to partially address the question of finding an optimal Mercer kernel for a given situation. The subsections~\ref{sec:specprob} and~\ref{sec:genform} discuss these theoretical aspects that are illustrated by a simple numerical simulation in Subsection~\ref{sec:Toy}.  A reader interested in the applications to dimensionality reduction can directly start the reading at Subsection~\ref{sec:Diff} where we anticipate the main results of our work which are applied to datasets in the Subsection~\ref{sec:Data}.
\item[-] In view of a specific choice of kernel related to a diffusion geometry and in order to make the paper more self-contained, we briefly review Diffusion Maps in Section~\ref{sec:DiffGeom}. Subsequently, in Section~\ref{sec:SDP}, the properties of a solution of the SDP that we propose to solve are discussed. The out-of-sample extension of the normalized diffusion kernel $K$  is explained in Section~\ref{sec:oos}.
\item[-]Finally, as it was mentioned hereinbefore, the numerical strategy and the simulations are reported in Section~\ref{sec:NumStrat}.
\end{itemize}

\section{From spectral problems to Semi-Definite Programs \label{sec:2}}
In this section, we first phrase our approach in a general measurable space in order to show that under certain hypotheses finding an optimal Mercer kernel makes sense.  This is in contrast with the more customary problem of learning a kernel matrix and may be called as ``learning the kernel''. Specifically, we consider integral operators with respect to a finite Borel measure until section~\ref{sec:Diff}. Then, we discuss the problem with respect to atomic measures supported on a finite set of points, so that the situation becomes finite dimensional and integrals become quadrature formulas.

Let us remind the classical setting of integral operators associated to continuous symmetric kernels.
For a thorough mathematical treatment, we refer to the book of Cucker and Zhou~\cite{CuckerZhou}. Formally, let $(X,\mu)$ be a measurable space where $\mu$ is a finite Borel measure (i.e., $\mu(X)<\infty$) and let $\mathcal{L}^2_\mu(X)$ be the space of square integrable functions.  Let $K:X\times X\to \mathbb{R}$ be a continuous symmetric kernel function. In view of applications, kernel functions are often taken to be continuous. Let us also assume \emph{a priori} that the measure $\mu$ is  supported on a subset $\supp\mu\subseteq X$. A natural object associated to $(X,\mu)$ and a kernel function is the integral operator $L_{K}:\mathcal{L}^2_\mu(X) \to \mathcal{C}(X)$ given by \[
(L_{K}f)(x) = \int_{X}K(x,y)f(y)\mathrm{d}\mu(y),\quad x \in X,
\]
where $f\in \mathcal{L}^2_\mu(X)$ is a real valued function on $X$. Since we have the injection $\mathcal{C}(X)\hookrightarrow \mathcal{L}^2_\mu(X)$, we think of $L_{K}$ as a linear operator valued in $\mathcal{L}^2_\mu(X)$. This integral operator is compact since its kernel is continuous (see Proposition 4.5, in~\cite{CuckerZhou}) and $L_{K}: \mathcal{L}^2_\mu(X)\to \mathcal{L}^2_\mu(X)$ is self-adjoint because $K$ is symmetric.
\begin{mydef}[Positive Semi-definite Kernel]
Let $\Omega\subseteq X$. The kernel $K$ is said to be positive semi-definite (p.s.d.) on $\Omega\times \Omega$ if, for all finite set $\{x_1, \dots ,x_m\}\subseteq \Omega$ with $m\in \mathbb{N}$, the kernel matrix whose entry $(i,j)$ is $ K(x_i,x_j)$ is positive semi-definite. Then, we write $K\succeq 0$.
\end{mydef}
In view of applications, the kernel $K$ is assumed to be \emph{p.s.d.} on $\Omega = \supp(\mu)$. This assumption is common in the context of kernel methods since such a kernel is related to a reproducing kernel Hilbert space.
\begin{mydef}[Mercer Kernel]
A continuous, symmetric and p.s.d. kernel function is called a Mercer kernel.
\end{mydef}Importantly, a Mercer kernel admits a countable orthonormal basis of continuous eigenfunctions $e_{\ell}\in \mathcal{L}^2_\mu(X)$ for $\ell \in \mathbb{N}$ with corresponding non-negative eigenvalues $\lambda_0\geq \lambda_1\geq \lambda_2\geq \dots$ sorted in decreasing order.
There is an absolutely convergent spectral decomposition $K(x,y) = \sum_{\ell\geq 0}\lambda_\ell e_\ell(x)e_\ell(y),$
of the Mercer kernel $K$. We have in particular the trace formula~\cite{Duflo,Brislawn} $\int_X K(x,x)\mathrm{d}\mu(x) = \sum_{\ell \geq 0}\lambda_\ell< +\infty$ and also, $\|L_{K} \|^2_{HS}=\sum_{\ell \geq 0}\lambda_\ell^2< +\infty$.
To sum up, $L_{K}$ is a Hilbert-Schmidt trace class operator. Let us use this eigenvector expansion in order to motivate our approach. Incidentally, notice also that any Hilbert-Schmidt operator can be represented as an integral operator.

\subsection{Spectral problems\label{sec:specprob}}
Major methods for dimensionality reduction such as Diffusion Maps or Kernel Principal Component Analysis rely on the calculation of the top eigenvectors of a kernel matrix. The eigenvector of maximal eigenvalue is denoted here by $e_0$. In the framework of this paper, the top eigenvector of $L_K$ is associated to the unnormalized projector
$
(L_{\rho^\star_0} f)(x) =\Tr(L_K)\times \int_{X}e_0(x)e_0(y)f(y)\mathrm{d}\mu(y)$, with $x \in X$, which  is an integral operator associated with the Mercer kernel $\rho^\star_0 = \Tr(L_K) e_0(x)e_0(y)$.
This finite rank operator is the solution of the problem for finding the supremum
\begin{equation}
\sup\Big\{ \big\langle L_{\rho}, L_{K}\big\rangle_{HS} : \rho\succeq 0 \text{ and } \Tr(L_\rho) =\Tr(L_K)\Big\}, \label{eq:FormalDefSDPTopEig}
\end{equation}
where the maximization is over Mercer kernels defined on $\Omega\times \Omega$ and where $\langle \cdot, \cdot\rangle_{HS}$ is the Hilbert-Schmidt inner product. Importantly, the fact that the supremum is attained is a non trivial question, in the infinite dimensional case  (We refer the reader to Proposition~\ref{Prop:HSop}).

Noticeably, the value of the  trace in the constraint $\Tr(\rho) =\Tr(K)$ is arbitrary since it only modifies the optimal solution by a constant normalization factor.
Then, the computation of the second top eigenvector of $L_{K}$ would require to solve another maximization problem involving the integral operator associated to the kernel $K^{(1)}(x,y) = K(x,y)-\lambda_0 e_0(x)e_0(y)$.
Clearly, by replacing $K$ by $K^{(1)}$ in~(\ref{eq:FormalDefSDPTopEig}), we obtain the following optimization problem for finding the supremum
\begin{equation}
\sup\Big\{ \big\langle L_{\rho},L_{K^{(1)}} \big\rangle_{HS} : \rho\succeq 0 \text{ and } \Tr(L_\rho) =\Tr(L_K)\Big\},  \label{eq:FormalDefSDPTopSecondEig}
\end{equation}
which yields the solution $
(L_{\rho^{(1)\star}} f)(x) =\Tr(L_K)\times e_1(x)\int_{X}e_1(y)f(y)\mathrm{d}\mu(y)$ with  $x \in X$.
Typically, a numerical method for calculating the top eigenvectors of a symmetric \emph{p.s.d.} matrix involves the power iteration algorithm. In the context of Diffusion Maps or spectral clustering, the quality of the spectral truncation strongly depends on the spectrum of the normalized kernel as illustrated for example in Figure~\ref{Fig:Favorable}. Indeed, a favorable case occurs when the spectrum contains several top eigenvalues close to $1$ followed by a large spectral gap. Clearly, spectral truncation is then a very reasonable approximation and one sees that only two eigenvectors are necessary to reduce the dimensionality of the dataset. Hence, this involves the computation of two top eigenvectors which amounts, e.g., to running twice the power iteration algorithm. However, it is intriguing to remark that the convergence rate of this algorithm typically depends on the ratio $\lambda_2/\lambda_1\approx 1$ for the calculation of $e_1$. Hence, this motivates the study of an alternative (non-spectral) problem for finding a potentially low rank \emph{p.s.d.} matrix incorporating the information contained in the subspace spanned only by the top eigenvectors.
\begin{figure}[h]
\setlength\figureheight{0.25\textwidth}
\setlength\figurewidth{0.25\textwidth}
\subfloat[][]{
\input{OriginalDataSpectrumPerfect.tikz}
}
\subfloat[][]{
\input{DMPerfect.tikz}}
\subfloat[][]{
\input{DMSpectrumPerfect.tikz}
}
\caption{Toy data with three clusters  (left). The embedding given by Diffusion Maps (middle) and the six largest eigenvalues of the kernel~(\ref{eq:KDiff}) with $\sigma =  1.5$ and, in particular, $\lambda_1 = 0.99$ and $\lambda_2 = 0.97$ (right). The use of only two eigenvectors for the embedding is clearly motivated by the spectrum. An illustration of the SDP embedding of a similar dataset can be found in Figure~\ref{Fig:OutliersBandwidth}.\label{Fig:Favorable}}
\end{figure}
\subsection{Towards Semi-definite Programs: general formulation \label{sec:genform}}
In order to avoid solving successive problems, let us now introduce another optimization problem which is the starting point for developing the novel contributions of this paper.

We define an approximation problem as follows. We want to find
\begin{equation}
s^\star = \sup \Big\{\big\langle L_{\rho}, L_{K}\big\rangle_{HS}: \rho\in \Gamma_K\Big\} ,\text{ with } \Gamma_K = \Big\{\rho \succeq 0: \rho(x,x) =K(x,x)\ {\rm for\ all}\ x\in \Omega\Big\}, \label{eq:FormalDefSDP}
\end{equation}
where the maximization is over all Mercer kernels with the same diagonal element as $K$. Remark that the diagonal of a \emph{p.s.d.} kernel is always non-negative. In addition, a convex combination of two \emph{p.s.d.} kernels with identical diagonal is still a \emph{p.s.d.} kernel with the same diagonal.
As a matter of fact, a feasible $L_\rho$ for~(\ref{eq:FormalDefSDP}) also satisfies the weaker condition $\|L_\rho\|_1 = \Tr(L_\rho) = \Tr(L_K)<\infty$, where $\|\cdot\|_1$ is the trace norm of the Banach space of trace-class operators. In other words, in view of the constraints, $L_\rho$ is constrained within a sphere in the space of trace class operators. In the same spirit, since $L_K$ is self-adjoint, the objective function is equivalently given by $\big\langle L_{\rho}, L_{K}\big\rangle_{HS} = \|L_{\rho}\circ L_{K}\|_1 = \Tr(L_{\rho}\circ L_{K})$. Let us list some reasons to propose the optimization problem~(\ref{eq:FormalDefSDP}).
\begin{itemize}
\item As a first motivation, the set of integral operators $\{L_\rho : \rho\in \Gamma_K\}$ is uniformly bounded since
\[
\|L_\rho\|_{op}\leq \|L_\rho\|_{HS}\leq \|L_\rho\|_{1}= \Tr(L_K),
\]
for all $\rho\in \Gamma_K$, where $\|L_\rho\|_{op}$ is the operator norm. The same is also true for the spectral approach.
\item A second motivation for proposing to maximize~(\ref{eq:FormalDefSDP}) over Mercer kernels with a fixed diagonal, is that all these kernel functions are uniformly bounded, i.e., for all $\rho\in \Gamma_K$ and for all $x$ and $y$ in $\Gamma_K$, we have
\begin{equation}
|\rho(x,y)|\leq \sqrt{K(x,x)} \sqrt{K(y,y)}\leq C_K^2, \  \text{ with } \ C_K = \sup_{x\in X}\sqrt{K(x,x)},\label{eq:UniBounded}
\end{equation}
as a consequence of the condition of positive definiteness. This is in contrast with the spectral approach.
\item Finally, if $K(x,x)$ is interpreted intuitively as the similarity of $x$ with itself, it is natural to preserve this property by requiring $\rho(x,x) = K(x,x)$ for all $x\in \Omega$.
\end{itemize}

As it is stated in Proposition~\ref{Prop:Bounded}, the objective function is trivially bounded on the set of kernels satisfying $\rho(x,x) = K(x,x)$ for all $x\in \Omega$. Therefore, the supremum of the objective in~(\ref{Prop:Bounded}) exists, i.e., is not infinite.
\begin{myprop}\label{Prop:Bounded}
Let $\rho$ and $K$ be Mercer kernels such that $\Tr(L_{\rho}) =\Tr(L_{K})$. Then, we have
\[
0\leq \big\langle L_{\rho}, L_{K}\big\rangle_{HS}\leq \lambda_0\Tr(L_{K})\leq \|L_{K} \|_{HS}\Tr(L_{K}),
\]
where $\lambda_0$ is the largest eigenvalue of $L_{K}$. In particular, by taking $\rho = K$, we also have the upper and lower bounds
$
\|L_K\|_{HS}^2\leq s^\star \leq \lambda_0\Tr(L_{K}),
$
where $s^\star$ is the supremum~(\ref{eq:FormalDefSDP}).
\end{myprop}
\begin{proof}
We have $\lambda_\ell(\rho)\leq \Tr(L_K)$ for all integers $\ell\geq 0$, and therefore,
\[
\|L_{\rho} \|^2_{HS} = \sum_{\ell\geq 0} \big(\lambda_\ell(\rho)\big)^2\leq\Tr(L_{K})\sum_{\ell\geq 0}\lambda_\ell(\rho) =  \big[\Tr(L_{K})\big]^2.\] Then, by the Cauchy-Schwarz inequality, it holds $
0\leq \big\langle L_{\rho}, L_{K}\big\rangle_{HS}\leq \|L_{K} \|_{HS}\|L_{\rho} \|_{HS}\leq\|L_{K} \|_{HS}\Tr(L_{K}).$
\end{proof}
In order to better understand the problem~(\ref{eq:FormalDefSDP}), we provide two simple instances when the supremum is reached.
A very elementary example saturating the first bound of Proposition~\ref{Prop:Bounded}, such that the supremum is attained, is given in Remark~\ref{Rem:1}. Next, we explain in Remark~\ref{Rem:Positivity} how the supremum is attained when the initial kernel $K$ takes only positive values.
\begin{myremark}[Tightness for the rank one case]\label{Rem:1}
Let $\mu$ be the Lebesgue measure on the compact subset $\Omega\subset \mathbb{R}^d$. Let $K$ the Mercer kernel $K(x,y) = \varphi(x)\varphi(y)$ for all $x,y\in \Omega$, where $\varphi$ is a continuous function $\varphi\in \mathcal{L}_\mu^2(\Omega)$ which satisfies $\Tr(L_K) = \int_\Omega \varphi(x)^2\mathrm{d}\mu(x)<\infty$. This choice includes as a particular case the linear kernel used in machine learning tasks. Trivially, the unique non-zero eigenvalue of K is $\lambda_0 = \Tr(L_K)$. Next, we simply find that the supremum of~(\ref{eq:FormalDefSDP}) is attained at $\rho^\star(x,y) = K(x,y)$, Indeed, we calculate
\[
\big\langle L_{\rho^\star}, L_{K}\big\rangle_{HS} = \int_\Omega K(x,y)\rho^\star(x,y)\mathrm{d}\mu(x)\mathrm{d}\mu(y) = \int_\Omega \varphi(x)^2\mathrm{d}\mu(x)\int_\Omega \varphi(y)^2\mathrm{d}\mu(y) = \lambda_0 \Tr(L_K),
\]
which is indeed the best upper bound as it is given in Proposition~\ref{Prop:Bounded}. This trivial case is naturally identical to the spectral approach.
\end{myremark}
\begin{myremark}[Exact solution for $K(x,y)\geq 0$]\label{Rem:Positivity}
In view of the theory of Hilbert-Schmidt operators, the objective function involves the integral
\[
\big\langle L_{\rho}, L_{K}\big\rangle_{HS}= \int_\Omega K(x,y)\rho(x,y)\mathrm{d}\mu(x)\mathrm{d}\mu(y)\geq 0.
\]
Since $\rho(x,y)$ is a \emph{p.s.d.} kernel, then for all $x$ and $y$ in $X$, we have $|\rho(x,y)|^2\leq  \rho(x,x)\rho(y,y)$ (see also~(\ref{eq:UniBounded})). In general,
the objective function is bounded as follows
\[
\big\langle L_{\rho}, L_{K}\big\rangle_{HS}\leq\int_\Omega | K(x,y)\rho(x,y)|\mathrm{d}\mu(x)\mathrm{d}\mu(y)\leq \int_\Omega |K(x,y)|\sqrt{K(x,x)}\sqrt{K(y,y)}\mathrm{d}\mu(x)\mathrm{d}\mu(y),
\]
for all $\rho\in \Gamma_K$. As a consequence, if we suppose that $K(x,y)\geq 0$ for all $x$ and $y$ in $X$,  the maximum of the objective is attained for the Mercer kernel \[\rho^\star (x,y) = \sqrt{K(x,x)}\sqrt{K(y,y)},\]
which is an explicit solution, although it is of little practical interest. This remark motivates why we only deal with kernels $K$ taking both positive and negative values in this paper.
\end{myremark}

We now state a result about the existence of a Hilbert-Schmidt operator attaining the supremum~(\ref{eq:FormalDefSDP}). For simplicity, we introduce several notations. Let $S(\mathcal{L}_\mu^2(X))$ be the Hilbert space of Hilbert-Schmidt operators on $\mathcal{L}_\mu^2(X)$. We also introduce the ball $B = \{L \in S(\mathcal{L}_\mu^2(X)): \|L\|_{HS}\leq \Tr(L_K)\}$ in the space of Hilbert-Schmidt operators. Let us first give elementary results.
\begin{mylemma}\label{Lemma:HS}
Let $g_1\in \mathcal{L}_\mu^2(X)$ and $g_2\in \mathcal{L}_\mu^2(X)$. Then, the integral operator $T_{g_1 g_2}:\mathcal{L}_\mu^2(X)\to \mathcal{L}_\mu^2(X)$ given by
\[
(T_{g_1 g_2}f)(x) = \int_{X}g_1(x)g_2(y)f(y)\mathrm{d}\mu(y),\qquad x\in X,
\]
is a Hilbert-Schmidt operator. Moreover, for all Hilbert-Schmidt $T\in S(\mathcal{L}_\mu^2(X))$, we have the identity $\langle  T_{g_1 g_2},T\rangle_{HS} = \langle g_1, Tg_2\rangle_{\mathcal{L}_\mu^2}$.
\end{mylemma}
\begin{proof}
Let $\{e_\ell\}_{\ell\geq 0}$ be an  orthonormal basis  of $\mathcal{L}_\mu^2(X)$.
We simply calculate the Hilbert-Schmidt norm as follows
\[
\|T_{g_1 g_2}\|^2_{HS} = \sum_{\ell\geq 0}\langle T_{g_1 g_2} e_\ell,T_{g_1 g_2} e_\ell\rangle_{\mathcal{L}_\mu^2} =  \langle g_1 ,g_1\rangle_{\mathcal{L}_\mu^2}\sum_{\ell\geq 0}\langle g_2 ,e_\ell\rangle_{\mathcal{L}_\mu^2}\langle e_\ell ,g_2\rangle_{\mathcal{L}_\mu^2} = \langle g_1 ,g_1\rangle_{\mathcal{L}_\mu^2}\langle g_2 ,g_2\rangle_{\mathcal{L}_\mu^2}<+\infty,
\]
which shows that $T_{g_1 g_2}$ is Hilbert-Schmidt. Then, the identity $\langle  T_{g_1 g_2},T\rangle_{HS} = \langle g_1, Tg_2\rangle_{\mathcal{L}_\mu^2}$ for all $T\in S(\mathcal{L}_\mu^2(X))$ follows from the definition of the Hilbert-Schmidt inner product.
\end{proof}

\begin{myprop}[Optimal HS operator]\label{Prop:HSop}
Let $\Gamma_K$ be the set of Mercer kernels as it is defined in~(\ref{eq:FormalDefSDP}), and let $\{\rho_n\}_{n\in \mathbb{N}}$ be a sequence of Mercer kernels in $\Gamma_K$ such that $\lim_{n\to +\infty}\big\langle L_{K}, L_{\rho_{n}}\big\rangle_{HS} = s^\star$, where $L_{\rho_{n}}\in B$ for all $n\in \mathbb{N}$.
Then, there exists a subsequence $\{\rho_{n_k}\}_{k\in \mathbb{N}}$ and a Hilbert-Schmidt operator $L^\star\in B$ such that
\[
\lim_{k\to+\infty}\big\langle L,L_{\rho_{n_k}} \big\rangle_{HS} = \big\langle L,L^\star \big\rangle_{HS}
, \text{ for all } L\in S(\mathcal{L}_\mu^2(X)).\]
In particular, by taking $L = L_K$ it holds that $
s^\star =\lim_{k\to+\infty}\big\langle L_K,L_{\rho_{n_k}} \big\rangle_{HS} = \big\langle L_K, L^\star\big\rangle_{HS}$.
\end{myprop}
\begin{proof}
We know that there is a sequence of $\rho_n$ such that $\big\langle L_{K},L_{\rho_{n}} \big\rangle_{HS} \to s^\star$ as $n\to \infty$, where each $L_{\rho_n}$ is in the ball $B$ in the Hilbert space $S(\mathcal{L}_{\mu}^2)$.
By the Banach-Alaoglu theorem, a closed ball in a Hilbert space is weakly compact, i.e., all sequence in the ball admits a subsequence which converges weakly. This property applied to the sequence $\{L_{\rho_{n}}\}_{n\in \mathbb{N}}$ yields the desired result.
\end{proof}
The Hilbert-Schmidt operator obtained thanks to the weak compactness is associated to a kernel function with good properties given in Corollary~\ref{Coral:Optrhostar}.
\begin{mycorol}[Optimal p.s.d symmetric kernel]~\label{Coral:Optrhostar}
There exists an integral kernel $\rho^\star\in \mathcal{L}_{\mu\otimes \mu}^2(X \times X)$ such that
\[
(L^\star f)(x) = \int_{X}\rho^\star(x,y)f(y)\mathrm{d} \mu(y), \ x\in X,
\]
and $\|L^\star\|_{HS} = \int_{X}|\rho^{\star}(x,y)|^2\mathrm{d} \mu(x)\mathrm{d} \mu(y)<\infty$.  Furthermore, for all functions $f$ and $g$ in $\mathcal{L}_\mu^2(X)$, we have
\begin{enumerate}
\item[(i)] $\langle f, L^\star f\rangle_{\mathcal{L}_\mu^2} =\int_{X}\rho^{\star}(x,y)f(x)f(y)\mathrm{d} \mu(x)\mathrm{d} \mu(y) \geq 0$,
\item[(ii)]$\langle f, L^\star g\rangle_{\mathcal{L}_\mu^2} = \langle  L^\star f, g\rangle_{\mathcal{L}_\mu^2},$
\item[(iii)]$\Tr(| L^\star|) = \Tr(L^\star)\leq \Tr(L_K)$
\end{enumerate}
so that $\rho^\star(x,y)$ is symmetric and positive semi-definite.
\end{mycorol}
\begin{proof}
We first recall that finite rank operators are dense in Hilbert-Schmidt operators. Hence, any Hilbert-Schmidt operator $T\in S(\mathcal{L}_\mu^2(X))$ is an integral operator with a kernel in $k\in \mathcal{L}_{\mu\otimes \mu}^2(X \times X)$. Indeed, let $\{e_\ell\}_{\ell\geq 0}$ be an orthonormal basis of $\mathcal{L}_\mu^2(X)$, then the kernel is given by
\[
k(x,y) = \sum_{\ell\geq 0} (T e_\ell)(x) e_\ell(y),
\]
as a result of Hilbert-Schmidt operator theory.  Then, the property (i) is obtained by choosing for $L$ in Proposition~\ref{Prop:HSop} the integral operator $T_{f^2}$ of kernel $f(x)f(y)$ thanks to Lemma~\ref{Lemma:HS}. Indeed, for all $k\in \mathbb{N}$,
$
\big\langle L_{\rho_{n_k}}, L_{f^2}\big\rangle_{HS}\geq 0
$ and therefore, by taking the limit, $\big\langle L^\star, L_{f^2}\big\rangle_{HS}\geq 0$. Still by Lemma~\ref{Lemma:HS}, $L^\star$ is a positive semi-definite integral operator: $\langle f, L^\star f\rangle_{\mathcal{L}_\mu^2}\geq 0$. Similarly, by choosing $L =T_{fg}$ and $L=T_{gf}$ in Proposition~\ref{Prop:HSop}  as defined in Lemma~\ref{Lemma:HS}, we find for each term in the subsequence of Proposition~\ref{Prop:HSop}
\[
\big\langle L_{fg} , L_{\rho_{n_k}}\big\rangle_{HS} = \big\langle L_{gf}, L_{\rho_{n_k}}\big\rangle_{HS},
\]
since each $L_{\rho_{n_k}}$ is a symmetric operator.
Then, by taking the limit on both sides and since the functions are all real valued, we find \[\langle f, L^\star g\rangle_{\mathcal{L}_\mu^2} = \big\langle L_{fg}, L^\star\big\rangle_{HS} = \big\langle L_{gf},L^\star \big\rangle_{HS}=\langle g, L^\star f\rangle_{\mathcal{L}_\mu^2}=\langle  L^\star f,g\rangle_{\mathcal{L}_\mu^2},\]
and therefore $L^\star$ is a symmetric operator and $\rho^\star(x,y)$ is a symmetric kernel.
In order to prove (iii), we select an orthonormal basis $\{e_\ell\}_{\ell\geq 0}$  of $\mathcal{L}_\mu^2(X)$. Then, by the same argument as before, we have
\[
\lim_{k\to \infty}\langle e_\ell, L_{\rho_{n_k}}e_\ell\rangle_{\mathcal{L}_\mu^2} = \langle e_\ell, L^\star e_\ell\rangle_{\mathcal{L}_\mu^2}=\lim_{k\to \infty}\inf \langle e_\ell, L_{\rho_{n_k}}e_\ell\rangle_{\mathcal{L}_\mu^2},
\]
for all $\ell\in \mathbb{N}$. Furthermore, notice that $\langle e_\ell, L_{\rho_{n_k}}e_\ell\rangle_{\mathcal{L}_\mu^2}\geq 0$ for all $\ell$ and $k\in\mathbb{N}$, and that
\[
\lim_{k\to \infty}\sum_{\ell\in \mathbb{N}}\langle e_\ell, L_{\rho_{n_k}}e_\ell\rangle_{\mathcal{L}_\mu^2} =\Tr(L_{K}) = \lim_{k\to \infty}\inf\sum_{\ell\in \mathbb{N}}\langle e_\ell, L_{\rho_{n_k}}e_\ell\rangle_{\mathcal{L}_\mu^2}.
\]
In consequence, by using Fatou's Lemma, we obtain
\[
\Tr(L_{K}) = \lim_{k\to \infty}\inf\sum_{\ell\in \mathbb{N}}\langle e_\ell, L_{\rho_{n_k}}e_\ell\rangle_{\mathcal{L}_\mu^2}\geq \sum_{\ell\in \mathbb{N}}\lim_{k\to \infty}\inf\langle e_\ell, L_{\rho_{n_k}}e_\ell\rangle_{\mathcal{L}_\mu^2} = \sum_{\ell\in \mathbb{N}}\langle e_\ell, L^\star e_\ell\rangle_{\mathcal{L}_\mu^2}.
\]
We now recall that $L^\star$ is a positive semi-definite and symmetric integral operator. Therefore, it holds that $|L^\star| = L^\star$ and $\Tr(L^\star)\leq \Tr(L_K)$.
\end{proof}
Let us choose the Lebesgue measure for $\mu$ and $X \subset \mathbb{R}^d$. A first natural question related to Corollary~\ref{Coral:Optrhostar} concerns the diagonal elements of the integral kernel of $L^\star$. Hence, strictly speaking, although the kernels in $\Gamma_K$ have a well-defined diagonal, the diagonal elements of the integral kernel $\rho^\star (x,y)$ constitute sets of Lebesgue measure zero. However, by using similar arguments than in the proof of Corollary~\ref{Coral:Optrhostar}, we can show that the average value of the integral kernel $\rho^\star (x,y)$ in a ball centered about a diagonal element is preserved.
Let $\rho\in \Gamma_K$ be a kernel in the converging subsequence of Proposition~\ref{Prop:HSop} and $x_0\in X$. Let $\epsilon>0$ small enough such that the ball of radius $\epsilon$ centered at $x_0\in X$  is in $X$, that is, $B_\epsilon(x_0)\subseteq X$. We denote by $1_{B_\epsilon(x_0)}(\cdot)$ the indicator function of this ball. Then, we have
\[
\int_{X}\rho^{\star}(x,y)1_{B_\epsilon(x_0)}(x)1_{B_\epsilon(x_0)}(y)\mathrm{d} x\mathrm{d} y = \int_{X}\rho(x,y)1_{B_\epsilon(x_0)}(x)1_{B_\epsilon(x_0)}(y)\mathrm{d} x\mathrm{d} y.
\]
A second issue is that the integral kernel $\rho^\star (x,y)\in\mathcal{L}_{\mu\otimes \mu}^2(X \times X)$ is not necessarily equivalent to a continuous function on $X\times X$. However, another notion of ``continuity'' of integral kernels associated to nuclear operators has been studied in~\cite{Vershik2013,Vershik2015}. Since the associated integral operator $L^\star$ is nuclear (or trace class), it can be shown that the integral kernel $\rho^\star (x,y)$ is ``virtually continuous''. The definition of virtually continuous measurable functions of two variables given in~\cite{Vershik2015} generalizes the ``almost continuous'' functions of one variable appearing in a classical theorem by Luzin.
In order to visualize better a numerical approximation of $\rho^\star (x,y)$ in a simple case, an illustration of a numerical solution of the optimization problem~(\ref{eq:FormalDefSDP}) is given in Section~\ref{sec:Toy}.

By defining an optimization problem over a more restrictive set of Mercer kernels, it may be possible to prove a uniform convergence to an optimal kernel as it is outlined in Remark~\ref{Rem:equicontinuity}.
\begin{myremark}\label{Rem:equicontinuity}
An alternative optimization problem can be phrased for instance by requiring that the Mercer kernels $\rho$  in $\Gamma_K$ are all Lipschitz continuous with the same Lipschitz constant (i.e., equicontinuous).
Then, since the set $\Gamma_K$ is uniformly bounded by~(\ref{eq:UniBounded}), we may apply the Arzela-Ascoli theorem (see, e.g.,~\cite{rudin}) in order to prove that the supremum is attained by a kernel in $\Gamma_K$.
\end{myremark}
In view of applications and numerical simulations, we now restrict the discussion to a finite dimensional Hilbert space.

\subsection{Towards Semi-definite Programs: formulation in finite dimension \label{sec:Diff}}
In the context of this paper and in view of diffusion geometry, we can suppose that $\Omega=\{x_1,\dots, x_N\}\subset X\subseteq\mathbb{R}^d$ is discrete.
Let us now assume that the measure $\mu$ is atomic and supported on the finite set of data points \[\Omega = \supp\mu=\{x_1,\dots, x_N\}\subset X.\]
We make a series of hypotheses on the kernel $K$.
\begin{assumption}Let $\supp\mu = \Omega\subseteq X$. We assume:
\begin{itemize}
\item[-] $K|_{\Omega\times\Omega}:\Omega\times \Omega\to \mathbb{R}$ is a Mercer kernel,
\item[-] for all $x\in X$, we have $K(x,x)>0$.
\end{itemize}\label{Assumption1}
\end{assumption}
Then, a specific choice of kernel $K$ satisfying Assumption~\ref{Assumption1} is done, i.e.,
\begin{equation}
K(x,y) = \frac{k(x,y)}{\sqrt{\bar{d}_\Omega(x)}\sqrt{\bar{d}_\Omega(y)}}-\sqrt{\frac{\bar{d}_\Omega(x)}{\vol(\Omega)}}\sqrt{\frac{\bar{d}_\Omega(y)}{\vol(\Omega)}}\label{eq:KDiff}
\end{equation}
for all $x,y\in X$ where the key hypothesis is that $k$ is a Mercer kernel taking \emph{strictly positive} values $k(x,y)> 0$ for all $x,y\in X$, while $\bar{d}_\Omega(x) =\sum_{z\in \Omega}k(x,z)$ and $\vol(\Omega) = \sum_{x\in \Omega}\bar{d}_\Omega(x)$. Notice that $\vol(\Omega)$ is not the volume in the sense of measure theory. For additional details and proofs, we refer to Section~\ref{sec:DiffGeom} and Section~\ref{sec:oos}.

Clearly, the formulation~(\ref{eq:FormalDefSDP}) bears a strong resemblance to Semi-Definite Programs used in the context of graph partitioning~\cite{Burer2003}, community detection or synchronization problems~\cite{MontanariPNAS}. A noticeable difference with previous works which are graph based is that our setting deals with general data points in a measurable space and is formulated in the context of kernel methods. Furthermore, the constraints on the diagonal elements in~(\ref{eq:FormalDefSDP}) also differ from previous approaches since the diagonal is not fixed to a constant value in this paper. These differences being clarified, we anticipate that, from a practical perspective, the numerical approach to solve~(\ref{eq:FormalDefSDP}) is not novel and based on a projected power method adapted from~\cite{BoumalGeNPowMeth}. Let us emphasize that the SDP is not solved directly but rather by following a method inspired by Burer and Monteiro~\cite{Burer2003}. Indeed, interior point methods for solving SDP may be difficult to scale up to large matrix sizes.

In this discrete context, the problem~(\ref{eq:FormalDefSDP}) reduces to the semi-definite program
\begin{equation}
\max_{\rho \succeq 0}\Tr (\rho K) ,\ {\rm subject\ to}\ \rho(x,x) =K(x,x)\ {\rm for\ all}\ x\in \Omega, \label{eq:OptimalityProperty}
\end{equation}
where the solution is a \emph{p.s.d.} kernel $\rho^\star:\Omega\times \Omega\to \mathbb{R}$. As a matter of fact, because of the choice of an atomic measure, the maximization problem is now over the set of all \emph{p.s.d.} symmetric matrices of fixed diagonal elements which constitute a convex set closely related to the so-called elliptope~\cite{LowRankSDP} (or the set of correlation matrices). Also, notice that, in principle, the discrete problem~(\ref{eq:OptimalityProperty}) makes sense for any continuous symmetric kernel with a strictly positive diagonal. Therefore, the constraint of positive definiteness of $K$ may also be relaxed, however we do not investigate this possibility here.
It is worth emphasizing that the rank of a solution $\rho^\star$ of the problem~(\ref{eq:OptimalityProperty}) is unconstrained but will be smaller than $N$. Indeed, from the number of constraints in the SDP~(\ref{eq:OptimalityProperty}), we already know (see, for instance, the work of Pataki~\cite{Pataki} or Theorem 1 in~\cite{Burer2003}) that the rank $r$ of a solution satisfies $
r(r+1)/2\leq N$.
Nevertheless, it is empirically observed in our simulations and in the literature~\cite{Burer2003,MontanariPNAS,LowRankSDP} that the solution of~(\ref{eq:OptimalityProperty}) has a rank much smaller than the previous bound. The low rank property of this solution may be understood from Proposition~\ref{Prop:Rank} that we prove in the sequel.

We anticipate now the main results of this paper and postpone the proofs to the subsequent sections.
Let $\ddiag(M)$ be a matrix of the same size as $M$ with zeros in its off-diagonal entries and the diagonal elements of $M$.
As shown hereafter, solving~(\ref{eq:OptimalityProperty}) yields an optimal kernel $ \rho^\star$ satisfying two conditions phrased in terms of the Laplacian-like matrix\begin{equation}
L(\rho^\star) = \ddiag(K)^{-1}\ddiag(K\rho^\star)-K.\label{eq:DualCertificate}
\end{equation}
Namely, the duality theory for semi-definite programs yields two conditions:
\begin{align}
&L(\rho^\star)\rho^\star = 0, &\text{ (complementary slackness) }\\ \label{eq:CompSlackness}
&L(\rho^\star)\succeq  0,  &\text{ (related to the dual feasibility)}
\end{align}
which are together sufficient for checking the optimality of a feasible $\rho^\star$ (cfr. Proposition~\ref{Prop:DualCertif}).
These conditions will prove themselves to be decisive in order to define an out-of-sample extension formula.

Intuitively, the combinatorial Laplacian has same structure as the certificate~(\ref{eq:DualCertificate}). Indeed, given a symmetric weight matrix $W$ with non-negative entries, the combinatorial Laplacian is the difference $D-W$, where  $D$ is the diagonal ``degree'' matrix with elements $D_{ii} = \sum_{j}W_{ij}$. For a more precise explanation of this analogy, we refer to Lemma~\ref{Lem:Trivial} and its proof.

The properties satisfied by the SDP embedding are the following:
\begin{itemize}
\item {\bf Rigidity}: a spectral decomposition of the solution of~(\ref{eq:OptimalityProperty}) yields $\rho^\star(x,y) = \sum_{\ell =1}^r\chi_\ell(x)\chi_\ell(y)$ for all $x,y\in \Omega$, where  $\chi_\ell$ is the eigenvector of $\rho^\star$ of eigenvalue $\lambda_\ell$ normalized such that $\sum_{x\in \Omega}\chi_\ell^2(x) = \lambda_\ell$. Hence, the mapping $\Xi:\Omega \to \mathbb{R}^r$ given by SDP embedding
\begin{equation}
x\mapsto \Xi(x) = \Big(
\chi_1(x),\dots,\chi_r(x)\Big)^\intercal,\label{eq:ApproxDiffusionMaps}
\end{equation}
which is  ``rigid'' since the length of the embedding vector is fixed as a constraint: $\|\Xi(x) \|_2  = \sqrt{K(x,x)}$ for all $x\in \Omega$. Hence, this guarantees that the data set will be embedded within a spherical shell
$\mathcal{S}_\Omega = \{y\in \mathbb{R}^r\ {\rm s.t.}\ \min_{z\in \Omega}\sqrt{K(z,z)}\leq \|y\|_2\leq \max_{z\in \Omega}\sqrt{K(z,z)}\}$.
\item {\bf Mean value property}: let $x\in \Omega$, then the identity
\begin{equation}
\chi_\ell(x) = \frac{K(x,x)}{\sum_{z\in \Omega}K(x,z)\rho^\star(z,x)}\sum_{y\in \Omega}K(x,y)\chi_\ell(y),\label{eq:Out-of-sample}
\end{equation}
holds as a consequence of the complementary slackness condition~(\ref{eq:CompSlackness}) for $\rho^\star$.   This property is indeed analogous to a mean value property since the positive multiplicative factor of the rhs of~(\ref{eq:Out-of-sample}) is smaller than $1/K(x,x)$ as stated in Lemma~\ref{Lem:Pos}. By analogy with the geometric harmonics~\cite{GeomHarm}, the identity~(\ref{eq:Out-of-sample}) can be used in order to define an extension of $\rho^\star$ on $X\subseteq \mathbb{R}^d$ as stated in~Theorem~\ref{Thm1} given hereafter. Another consequence of~(\ref{eq:CompSlackness}) is that the multiplication by $K$ defined by
\[
K\Xi(x) = \Big(
\sum_{y\in \Omega}K(x,y)\chi_1(y),\dots,\sum_{y\in \Omega}K(x,y)\chi_r(y)\Big)^\intercal,
\]
is ``conformal'', i.e., it preserves the angles between every pair of embedding vectors $\Xi(x_i)$ and $\Xi(x_j)$ for all $i$ and $j\in\{1,\dots, N\}$.

\item {\bf Out-of-sample extension} (projected Nystr\"om extension): Let $\ell\in\{1,\dots,r\}$.
The function $\bar{\chi}_\ell:X\to \mathbb{R}$ given by
\begin{equation}
\bar{\chi}_\ell(x) = \sqrt{\frac{K(x,x)}{\sum_{\ell=1}^{r}\Big(\sum_{y\in \Omega}K(x,y)\chi_\ell(y)\Big)^2}}\sum_{y\in \Omega}K(x,y)\chi_\ell(y), \text{ for all } x\in X,\label{eq:Out-of-sample2}
\end{equation}
is a bona fide extension of $\chi_\ell: \Omega\to \mathbb{R}$ as stated in Theorem~\ref{Thm:OOS}, i.e., the extension  $\bar{\chi}_\ell$ reduces to $\chi_\ell$ on $\Omega$.
 Notice that by construction the normalization ensures that $\sum_{\ell = 1}^{r}\bar{\chi}^2_\ell(x) = K(x,x)$, for all $x\in X$. In the sequel, the extension of a function will always be denoted by using an overbar.
\end{itemize}
\subsection{Out-of-sample extension of the embedding coordinates}
Concerning the out-of-sample extension, the scheme $\bar{\chi}_\ell(x)  = N(x) \sum_{y\in \Omega}K(x,y)\chi_\ell(y)$  with the constant factor $N(x) = 1/\lambda_\ell$ is typical of the  ``Nystr\"om approximation'' and defines the so-called geometric harmonics.
Let us also recall that any eigenvector $\psi_\ell:\Omega\to \mathbb{R}$ of $K$ can be extended to $\bar{\psi_\ell}:X\to \mathbb{R}$ thanks to a Nystr\"om extension formula~\cite{WilliamSeeger}
\[
\bar{\psi_\ell}(x) = \frac{1}{\lambda_\ell}\sum_{z\in \Omega} K(x,z)\psi_\ell(z), \quad x\in X,
\]
where the multiplicative factor is simply the inverse of the eigenvalue associated to $\chi_\ell$. Notice that the function given by the rhs of~(\ref{eq:Out-of-sample2}) can always be defined if $\chi_\ell$ is an eigenvector of any \emph{p.s.d.} matrix. However, given a positive semi-definite $\rho$, the formula~(\ref{eq:Out-of-sample2}) does not in general yield an extension of the eigenvector $\chi_\ell$. The difference with respect to any such kernel is that the solution $\rho^\star$ satisfies an optimality condition of the problem~(\ref{eq:OptimalityProperty}). This optimality condition happens to be crucial in order to define a bona fide out-of-sample extension, i.e., an extension which reduces to the original function when restricted to $\Omega$. Furthermore, contrary to the ``geometric harmonics'' method, the extension formula~(\ref{eq:Out-of-sample}) is non-linear in terms of the eigenfunction $\chi_\ell$. Hence, we can call the extension formula~(\ref{eq:Out-of-sample}) a ``projected Nystr\"om extension''.
\begin{figure}
\centering
\begin{tikzpicture}
  \shade[ball color = gray!40, opacity = 0.1] (0,0) circle (2cm);
  \draw (0,0) circle (2cm);
  \draw (-2,0) arc (180:360:2 and 0.6);
  \draw[dashed] (2,0) arc (0:180:2 and 0.6);
  \fill[fill=black] (0,0) circle (1pt);
  \draw[dashed] (0,0 ) -- node[below]{\small $K(\bar{x},\bar{x})$} (2,0);
    \draw[thick] (-1,1) -- (-2.5,2.5);
    \draw[dashed] (0,0) -- (-1,1);
    \draw (-2.5,2.5) node[left]{$\sum_{i = 1}^NK(\bar{x},x_i)\Xi(x_i)$};
    \draw (-2.5,3) node[left]{Nystr\"om extension};

\draw (-2.5,2.5) node{\textbullet};
 \draw (-1.,1.) node{\textbullet};
       \draw (-1.+0.1,1) node[right]{$\bar{\Xi}(\bar{x})$};
\draw (-1.+0.65,1.2) node[above]{projection};
 \draw[blue] (2,2) node{\textbullet};           \draw[blue] (2,2) node[right]{$\Xi(x_1)$};
 \draw[blue] (-2,-1.8) node{\textbullet};           \draw[blue] (-2.5,-1.8) node[above]{$\Xi(x_2)$};
 \draw[blue] (-0.8,-1.2) node{\textbullet};           \draw[blue] (-0.8+0.1,-1.2) node[right]{$\Xi(x_N)$};

\draw[blue] (-2.4,0.2) node{\textbullet};          \draw[blue] (-2.4,0.2) node[left]{$\Xi(x_3)$};
\draw[dotted,blue] (0,0) -- (0.83*-2.4,0.83*0.2);
\draw[blue](0.83*-2.4,0.83*0.2) -- (-2.4,0.2);

\draw[dotted,blue] (0,0) -- (-0.8,-1.2);
\draw[blue] (-2*0.6,-1.8*0.6) -- (-2,-1.8);
\draw[dotted,blue] (0,0) -- (-2*0.6,-1.8*0.6);

\draw[blue] (1,1) -- (2,2);
\draw[dotted,blue] (0,0) -- (1,1);
\end{tikzpicture}
\caption{Illustration of the projected Nystr\"om extension. The blue points represent the embedding of the initial data points in $\mathbb{R}^r$. The new point $\bar{x}$  in black is embedded by projecting on a sphere in $\mathbb{R}^r$ centered at the origin of radius $K(\bar{x},\bar{x})$ the point obtained thanks to the classical extension $\sum_{i = 1}^NK(\bar{x},x_i)\Xi(x_i)$. \label{Fig:Nystrom}}
\end{figure}
Let us now discuss the reasons why these properties hold.
First of all, an essential feature of the mean value property~(\ref{eq:Out-of-sample}) is that for all $x\in \Omega$ the proportionality factor $K(x,x)/\sum_{z\in \Omega}K(x,z)\rho^\star(z,x)$ is a real value between $0$ and $1/K(x,x)$.
\begin{mylemma}\label{Lem:Pos}
Let $x\in \Omega$. Then, it holds that $
\sum_{z\in \Omega}K(x,z)\rho^\star(z,x)\geq K^2(x,x)> 0$.
Hence, the diagonal matrix $\ddiag(K\rho^\star)$ is strictly positive definite.
\end{mylemma}
\begin{proof}
This is a consequence of the positive semi-definite property of $L(\rho^\star) = \ddiag(K)^{-1}\ddiag(K\rho^\star)-K$. Indeed, let $\delta_x$ be the unit delta function associated to $x\in \Omega$, then the result is simply given by the inequality $\delta_x^\intercal L(\rho^\star)\delta_x\geq 0$.
\end{proof}
It is now possible to prove the main theorem of this paper.  Indeed, Theorem~\ref{Thm:OOS} states that if $\chi_\ell:\Omega\to \mathbb{R}$ is an eigenvector of the optimal kernel $\rho^\star$ then the formula~(\ref{eq:Out-of-sample2}) gives rise to a well defined extension of the eigenfunction $\bar{\chi}_\ell:X \to \mathbb{R}$.
\begin{mythm}[Out-of-sample extension]\label{Thm:OOS}
Let $\chi_\ell$ for $\ell\in \{1,\dots,r\}$ be the eigenvectors of the solution $\rho^\star$ of~(\ref{eq:OptimalityProperty})  such that $\rho^\star(x,y) = \sum_{\ell=1}^{r}\chi_\ell(x)\chi_\ell(y)$. Then, the extension of the $\ell$-th eigenvector defined in~(\ref{eq:Out-of-sample2}) reduces to
$\bar{\chi}_\ell(x) = \chi_\ell(x)$,
for all $x\in \Omega$.
\end{mythm}
\begin{proof}
Firstly, by using~(\ref{eq:Out-of-sample}) and $\sum_\ell\chi_\ell^2(x) = K(x,x)$, we have
\begin{equation*}
\Big\{\sum_{\ell=1}^{r}\Big(\sum_{y\in \Omega}K(x,y)\chi_\ell(y)\Big)^2\Big\}^{-1/2} =\sqrt{ \frac{K(x,x)}{\Big(\sum_{z\in \Omega}K(x,z)\rho^\star(z,x)\Big)^2}}=\frac{\sqrt{K(x,x)}}{\Big|\sum_{z\in \Omega}K(x,z)\rho^\star(z,x)\Big|}.
\end{equation*}
Hence, by ~(\ref{eq:Out-of-sample}) and thanks to Lemma~\ref{Lem:Pos}, we find
\[
\bar{\chi}_\ell(x) =\sqrt{K(x,x)}\frac{\sqrt{K(x,x)}}{\Big|\sum_{z\in \Omega}K(x,z)\rho^\star(z,x)\Big|}\sum_{y\in \Omega}K(x,y)\chi_\ell(y) = \chi_\ell(x).
\]
\end{proof}
A natural question is: in what sense is the out-of-sample formula optimal? We partially address the question of the optimality of the extension of the embedding coordinates.

A major feature of the out-of-sample formula of the embedding coordinates is that it is optimal provided that the extension is considered only on one point. In order to formalize this result, let $X = \Omega \cup \{\bar{x}\}$ be the original space with one additional point. Then, following~(\ref{eq:Out-of-sample2}), we recall that  the value of the extended eigenfunction $\bar{\chi}_\ell: \Omega \cup \{\bar{x}\}\to \mathbb{R}$ is
\begin{equation}
\bar{\chi}_\ell(\bar{x}) =\sqrt{K(\bar{x},\bar{x})} \frac{\sum_{y\in \Omega}K(\bar{x},y)\chi_\ell(y)}{\sqrt{\sum_{\ell=1}^{r}\Big(\sum_{y\in \Omega}K(\bar{x},y)\chi_\ell(y)\Big)^2}}.\label{eq:ExtensionOnePoint}
\end{equation}
This definition yields an optimal value  of the extension in the following sense: it addresses the problem to maximize $\Tr(\bar{\rho}K)$ over all $\bar{\rho}:X\times X\to \mathbb{R}$ satisfying  $\bar{\rho}(x,y) = \sum_{\ell = 1}^r\chi_\ell(x)\chi_\ell(y)$ for all $x,y\in \Omega$ and $\bar{\rho}(\bar{x},\bar{x}) = K(\bar{x},\bar{x})$. Indeed, we propose an optimization problem which aims to find the best value of the extension given the eigenfunctions $\chi_\ell:\Omega\to \mathbb{R}$.
\begin{myprop}[Optimality of the extension]
The solution to
\[
\max_{u_\ell:X\to \mathbb{R}} \sum_{x,y\in \Omega\cup\{\bar{x}\}}K(x,y)\sum_{\ell =1}^{r}u_\ell(x)u_\ell(y) {\rm\ subject\ to } \begin{cases}
\sum_{\ell =1}^{r}u^2_\ell(\bar{x}) = K(\bar{x},\bar{x})\\
 u_\ell|_\Omega =\chi_\ell,
 \end{cases}
\] is
$u^\star_\ell(x) = \bar{\chi}_\ell(x)$ as given in~(\ref{eq:ExtensionOnePoint}).
\end{myprop}
\begin{proof}
We show that the optimal solution yields the desired extension.
Firstly, let us solve the constraints $u_\ell|_\Omega =\chi_\ell$. Then, the objective becomes
\[
\sum_{x,y\in \Omega}K(x,y)\sum_{\ell =1}^{r}\psi_\ell(x)\psi_\ell(y) +2\sum_{x\in \Omega}K(\bar{x},x)\sum_{\ell =1}^{r}\psi_\ell(x)u_\ell(\bar{x}) + K(\bar{x},\bar{x})\sum_{\ell =1}^{r}u^2_\ell(\bar{x}),
\]
where the first term is the optimal value of the objective in~(\ref{eq:OptimalityProperty}). Then, by the Lagrange multipliers technique, there exists $\Lambda\in \mathbb{R}$ such that for all $\ell = 1, \dots, r$, we have $\sum_{x\in \Omega}K(\bar{x},x)\psi_\ell(x)= \Lambda u_\ell(\bar{x})$
and $\sum_{\ell =1}^{r}u^2_\ell(\bar{x}) = K(\bar{x},\bar{x})$. We find
\[
\Lambda^{-2} = K(\bar{x},\bar{x})/\sum_{\ell =1}^{r}\big(\sum_{x\in \Omega}K(\bar{x},x)\psi_\ell(x) \big)^2,
\]
so that $\Lambda$ is determined up to a sign. The positive square root of $\Lambda^2$ maximizes the objective function and yields the desired out-of-sample formula.
\end{proof}

Actually, this optimality property is lost when the extension is done on several points, although the one-point extension can be used as an approximation.

\subsection{Out-of-sample extension of the kernel}
We firstly state the extension formula for the kernel defined on $\Omega\times\Omega$ and discuss hereafter the motivation of this choice.
The extension of the eigenvectors of the kernel $\rho^\star:\Omega\times\Omega\to \mathbb{R}$ with $\Omega\subseteq X$ yields an extension of the kernel itself to a \emph{p.s.d.} kernel $\bar{\rho}^\star:X\times X\to \mathbb{R}$.
\begin{mycorol}[Extended kernel via the extended eigenvectors]\label{Corol} The extension of the kernel $\rho^\star:\Omega\times \Omega\to \mathbb{R}$  given by the  kernel $\bar{\rho}^\star:X\times X\to \mathbb{R}$, such that \[\bar{\rho}^\star(x,x') =\sum_{\ell =1}^{r}\bar{\chi}_\ell(x)\bar{\chi}_\ell(x'),\]
for all $x,x'\in X$,  is positive semi-definite.
\end{mycorol}
\begin{proof}
This is a consequence of the definition of the formula~(\ref{eq:Out-of-sample2}).
\end{proof}

Hence, the extension of $\rho^\star$ is simply obtained thanks to the extension of its eigenvectors. For an explicit form of $\bar{\rho}$, we refer to~\ref{AppExtendedKernel}. Also, notice that the extended eigenvectors $\bar{\chi}_\ell$ are no more orthogonal with respect to each other.

We now argue that this strategy is essentially the most natural manner to extend $\rho^\star$. In order to answer this question, we consider the following  problem formulated in terms of matrices. Let us identify $\rho^\star$ with a $N\times N$ symmetric \emph{p.s.d.} matrix, while its basis of orthogonal eigenvectors $\chi_\ell$ are simply $N\times 1$ vectors normalized so that $\rho^\star = \sum_{\ell=1}^r \chi_\ell \chi_\ell^\intercal$, with $r\leq N$.
Then, the problem consists in finding the $(N+1)\times (N+1)$ matrix given by the block matrix
\begin{equation}
\bar{\rho}(b,s) =
\left[\begin{array}{c|c}
 \rho^\star & b\\
 \hline
 b^\intercal & s
\end{array}\right],\label{eq:mat}
\end{equation}
with $b\in \mathbb{R}^{N\times 1}$ and $s>0$. We now show that for any choice of $b\in \mathbb{R}^{N\times 1}$ in the range of $\rho^\star$ (with $\rho^{\star\intercal} =\rho^\star$), we can find $s>0$ such that
$\bar{\rho}(b,s)$ is \emph{p.s.d.}, whereas the minimal $s$ for which it is true yields an extension in terms of the extended eigenvectors.
We have Proposition~\ref{Prop:Extension}.
\begin{myprop}[Extended p.s.d. matrix]\label{Prop:Extension}
Let $\rho^\star\in \mathbb{R}^{N\times N}$ p.s.d. matrix with eigendecomposition $\rho^\star=\sum_{\ell=1}^{r} \chi_\ell \chi_\ell^\intercal$ and let $b\in \mathbb{R}^{N\times 1}$ be a vector such that $b\neq 0$. Then, the matrix $\bar{\rho}(b,s)\in \mathbb{R}^{(N+1)\times (N+1)}$ as given in~(\ref{eq:mat}) is \emph{p.s.d.} if and only if
\[b = \sum_{\ell=1}^{r} b_\ell \chi_\ell\in\range \rho^\star\ \text{ and }\ s\geq s_{\min} = \sum_{\ell = 1}^{r}b_\ell^2.\] In particular, for the lower bound $s = s_{\min}$, the extended matrix reads
$
\bar{\rho}(b,s_{\min}) = \sum_{\ell=1}^{r}\left[\begin{array}{c}
\chi_\ell \\
 \hline
  b_\ell
\end{array}\right]\left[\begin{array}{c |c}
 \chi_\ell^\intercal &  b_\ell
\end{array}\right].
$
\end{myprop}
\begin{proof}
Let
$\bar{v}^\intercal
= \left[
\begin{array}{c |c}
 v^\intercal &  v_s
\end{array}
\right]
\in\mathbb{R}^{1\times (N+1)}$, with $v\in \mathbb{R}^{N\times 1}$ and $v_s\in \mathbb{R}$. In this proof, we suppose that $v\neq 0$. In order to study the positivity of $\bar{\rho}$, we calculate $
O( v ,  v_s) \triangleq \bar{v}^\intercal\bar{\rho}(b,s)\bar{v}= v^\intercal\rho^\star v + 2v_s v^\intercal b + s v_s^2$ with $b\neq 0$. A lower bound on $O( v ,  v_s)$  is merely obtained by minimizing over all $v_s$ as follows:
\begin{equation}
O( v ,  v_s)\geq \min_{v_s\in \mathbb{R}}O( v ,  v_s) = O( v ,  v^\star_s) = v^\intercal(\rho^\star  - \frac{1}{s} b b^\intercal )v,\label{eq:LowerBound}
\end{equation}
with $v_s^\star= \arg\min_{v_s\in \mathbb{R}}O( v ,  v_s)= - b^\intercal v/s$.

Firstly, let us assume that $b\notin \range\rho^\star$. Then, we can write $b = b_0+ b'$ with $b_0 \in \ker \rho^\star$ and $b_0^\intercal b' =0$. Hence, by choosing $u = b_0\in \ker \rho^\star$ and by using~(\ref{eq:LowerBound}), we obtain
\[
O( u ,  u^\star_s) = u^\intercal(\rho^\star  - \frac{1}{s} b b^\intercal )u = -\frac{1}{s}(b_0^\intercal b_0)^2\leq 0,
\]
and therefore, $\bar{\rho}(b,s)$ is indefinite if $b\notin \range\rho^\star$.

Secondly, as a consequence, we assume that $b$ is in the range of $\rho^\star$, i.e.,  $b = \sum_{\ell =1}^{r} b_\ell \chi_\ell$. Let the components $v_\ell$ be defined such that $v = \sum_{\ell =1}^{r} v_\ell \chi_\ell/\lambda_\ell\neq 0$, where $\lambda_\ell$ is the eigenvalue associated with $\chi_\ell$. It holds that
\[
O( v ,  v^\star_s) = v^\intercal(\rho^\star  - \frac{1}{s} b b^\intercal )v = \sum_{\ell =1}^{r} v_\ell^2-\frac{1}{s}\Big(\sum_{\ell =1}^{r} b_\ell v_\ell\Big)^2.
\]
Then, if $s\geq s_{\min} = \sum_{\ell = 1}^{r}b_\ell^2$, we have the lower bound
\[
O( v ,  v^\star_s) = \sum_{\ell =1}^{r} v_\ell^2-\frac{1}{s}\Big(\sum_{\ell =1}^{r} b_\ell v_\ell\Big)^2\geq \Big(\sum_{\ell =1}^{r} v_\ell^2\Big)\Big[1-\frac{\big(\sum_{\ell =1}^{r}\hat{b}_\ell  v_\ell\big)^2}{\sum_{\ell =1}^{r} v_\ell^2}\Big], \text{ with } \hat{b}_\ell = \frac{b_\ell}{\sqrt{\sum_{\ell =1}^r b_\ell^2}}\leq 1.
\]
Since this is valid for all $v\in\mathbb{R}^{N\times 1}\setminus\{0\}$, we are going to study the worst case scenario, that is, when the lower bound reaches its minimum value. Indeed, we observe that the maximum of the Rayleigh quotient
\[
\max\Big\{ \frac{\big(\sum_{\ell =1}^{r}\hat{b}_\ell  v_\ell\big)^2}{\sum_{\ell =1}^{r} v_\ell^2}\Big| v_\ell\in\mathbb{R},\ \ell = 1,\dots, r \text{ and }\sum_{\ell =1}^{r} v_\ell^2\neq 0 \Big\} = 1,
\]
is obtained for $v_\ell^{\max} = \alpha \hat{b}_\ell$ for all $\alpha\neq 0$, which can be thought of as an eigenvector.
Therefore,
\[
O( v ,  v^\star_s)\geq\Big(\sum_{\ell =1}^{r} v_\ell^2\Big)\Big[1-\frac{\big(\sum_{\ell =1}^{r}\hat{b}_\ell  v_\ell\big)^2}{\sum_{\ell =1}^{r} v_\ell^2}\Big] \geq 0.
\]
We can now show that, if $s<\sum_{\ell = 1}^{r}b_\ell^2$ and $b\in \range \rho^\star$, then $\bar{\rho}$ is indefinite. Indeed, there exists $v^{(b)}$ such that  $O( v^{(b)} ,  v^{(b)\star}_s)<0$. In that case,  for all $v\neq 0$, we have
\[
O( v ,  v^\star_s) = \sum_{\ell =1}^{r} v_\ell^2-\frac{1}{s}\Big(\sum_{\ell =1}^{r} b_\ell v_\ell\Big)^2<\sum_{\ell =1}^{r} v_\ell^2-\Big(\sum_{\ell =1}^{r} \hat{b}_\ell v_\ell\Big)^2.
\]
By choosing $v = v^{(b)}$ such that $v^{(b)}_\ell = \hat{b}_\ell$, we find
\[
O( v^{(b)} ,  v^{(b)\star}_s) <\sum_{\ell =1}^{r} \hat{b}_\ell^2-\Big(\sum_{\ell =1}^{r} \hat{b}_\ell\hat{b}_\ell\Big)^2 = 0.
\]
Finally, in the special case $ s= s_{\min}=\sum_{\ell = 1}^{r}b_\ell^2$, we simply calculate that the extended \emph{p.s.d.} $\bar{\rho}(b,s_{\min})$ is simply obtained thanks to the extension of the eigenvectors $\bar{\chi}_\ell = \left[\begin{array}{c}
\chi_\ell \\
 \hline
  b_\ell
\end{array}\right]$ as follows \[
\bar{\rho}(b,s_{\min})=  \sum_{\ell=1}^{r}\left[\begin{array}{c|c}
 \chi_\ell \chi^\intercal_\ell&  b_\ell \chi_\ell\\
 \hline
  b_\ell \chi_\ell^\intercal &  b_\ell^2
\end{array}\right]= \sum_{\ell=1}^{r}\left[\begin{array}{c}
\chi_\ell \\
 \hline
  b_\ell
\end{array}\right]\left[\begin{array}{c |c}
 \chi_\ell^\intercal &  b_\ell
\end{array}\right].
\]
\end{proof}
Let us draw the attention on two facts. Firstly, the extension of the kernel is possible only if $b$ is restricted to be in the range of $\rho^\star$. This justifies why the extension formula in Corollary~\ref{Corol} includes only eigenvectors of non-zero eigenvalues. Secondly, in that case, the extended matrix with a minimal trace (or nuclear norm) is obtained thanks to the extension of the eigenvectors of $\rho^\star$.
Furthermore, if we define $\mathrm{k}_i = K(\bar{x}, x_i)$ for $i = 1, \dots, N$ and $\kappa = K(\bar{x},\bar{x})$, then, the extension proposed in~(\ref{eq:ExtensionOnePoint}) is given in the matrix notation~(\ref{eq:mat}) in terms of the kernel
\begin{equation}\label{eq:Kbar}
\bar{K} =
\left[\begin{array}{c|c}
 K & \mathrm{k}\\
 \hline
 \mathrm{k}^\intercal & \kappa
\end{array}\right],
\end{equation}
 as follows:
$
b = \sqrt{\frac{\kappa}{\mathrm{k}^\intercal \rho^\star \mathrm{k}}}\rho^\star \mathrm{k}$ or equivalently, $b_\ell = \sqrt{\frac{\kappa}{\mathrm{k}^\intercal \rho^\star \mathrm{k}}}\chi_\ell^\intercal \mathrm{k}$.
For a discussion of the relationship between $\bar{\rho}$ and the solution of an extended SDP in $N+1$ points, we refer the reader to~\ref{App:np1}.

In view of a practical implementation of our proposal, let us now quickly review the diffusion geometry setting which motivates the choice~(\ref{eq:KDiff}) for the kernel $K$.

\section{Diffusion geometry and Spectral truncation of the diffusion kernel\label{sec:DiffGeom}}
Let us restrict our discussion to the case where the data set is the discrete set $\Omega=\{x_1,\dots, x_N\}\subset X\subseteq\mathbb{R}^d$ endowed with a uniform atomic measure. The notion of similarity between two points is given here by a Mercer kernel $k$ taking only strictly positive values $k(x,y)> 0$ for all $x,y\in X$. An example is the Gaussian kernel
$k(x,y) = \exp(-\|x-y\|_2^2/\sigma^2)$ which is also a positive definite kernel, i.e., for any finite set $\Omega\subset \mathbb{R}^d$, its kernel matrix $[k(x,y)]_{x,y\in \Omega}$ is strictly positive definite.

In the context of diffusion geometry, it is customary to construct a reversible Markov chain on $\Omega$. The transition probability from $x$ to $y$ reads
\[
p(x,y) =  \frac{k(x,y)}{d_\Omega(x)},
\]
where the degree $d_\Omega(x) =\sum_{z\in \Omega}k(x,z)$ is a local measure of volume of $x\in \Omega$. Noticeably, the stationary distribution of the random walk at $x\in \Omega$ is proportional to the volume of the data point $x$, i.e.,
\[
\phi_0(x) = \frac{d_\Omega(x)}{\vol(\Omega)}>0,
\]
where the total volume is $\vol(\Omega) = \sum_{z\in \Omega}d_\Omega(z)$. Hence, if $k$ is the Gaussian kernel, the stationary distribution $\phi_0(x)$ is simply the Parzen window estimator of a probability density.

By construction, the random walk is reversible, i.e. $\phi_0(x)p(x,y) = \phi_0(y)p(y,x)$,  and  the eigenvalues of the transition matrix $\{\lambda_\ell\}_{\ell =0}^{N-1}$ are all positive, that is $1\geq \lambda_\ell\geq 0$ for $\ell =0,\dots, N-1$. Importantly, the positiveness of the eigenvalues is a consequence of the fact that $k$ is a \emph{p.s.d.} kernel. Incidentally, if $k$ is the Gaussian kernel, then $1\geq \lambda_\ell> 0$ for $\ell =0,\dots, N-1$. Indeed,  the transition matrix $p(x,y)$ is related by conjugation to the normalized kernel
\begin{equation}
k_\Omega(x,y) = \frac{k(x,y)}{\sqrt{d_\Omega(x)d_\Omega(y)}},\label{eq:NormalizedKernel}
\end{equation}
which is positive definite and shares the same eigenvector with the so-called normalized Laplacian is $L_N = \mathbb{I} -k_\Omega$. Notice that we use the subscript $\Omega$ to emphasize that the normalization is done with respect to the discrete set $\Omega$.
Since the random walk generated by $p(x,y)$ is reversible, the transition matrix at time $t\in \mathbb{N}_\star$ can be expanded in a bi-orthogonal basis as follows
\begin{equation}
p_{t}(x,y) =  \sum_{\ell =0}^{N-1}\lambda_{\ell}^t \psi_{\ell}(x)\phi_\ell(y),\label{eq:Transition Proba}
\end{equation}
with $\phi_\ell^\intercal\psi_{\ell'} = \delta_{\ell,\ell'}$, while $\phi_\ell(x) = \phi_0(x) \psi_\ell(x)$, for all $\ell\in \{0,\dots,N-1\}$. Then, the $\phi_\ell$'s are orthogonal with respect an inner product reweighted by the inverse volume $\phi^{-1}_0(y)$ of each data point $y$, that is,
\[
\langle f, g\rangle_{\phi^{-1}_{0}} = \sum_{y\in \Omega}f(y)g(y)\phi^{-1}_0(y),
\]
which also defines the corresponding norm given by $\|\cdot\|_{\phi_0^{-1}}$. Notice that~(\ref{eq:Transition Proba}) can also be defined for a continuous time $t>0$.
\begin{myremark}
The bi-orthogonal basis $(\phi_\ell,\psi_\ell)$ is related to an orthonormal basis $\{u_\ell\}_{\ell =0}^{N-1}$ by the relation $\phi_\ell(x) = \phi^{-1/2}_0(x) u_\ell(x)$ and $\psi_\ell(x) = \phi^{1/2}_0(x)u_\ell(x)$.
\end{myremark}

Still following~\cite{Coifman05geometricdiffusions}, the row $p_{t}(x,\cdot)$ is also interpreted as a bump function centered at $x$ of width depending on $t$. The $\ell_2$-distance between two bump functions at different data points, e.g. $p_{t}(x,\cdot)$ and $p_{t}(z,\cdot)$, is interpreted as the diffusion distance between the two data points .
More precisely, the diffusion distance is defined by the $\ell_2$-distance $D^2_t(x,z) = \|p_{t}(x,\cdot)-p_{t}(z,\cdot) \|^2_{\phi_0^{-1}}$.
In~\cite{Coifman05geometricdiffusions}, the diffusion map is defined as
\begin{align*}
\mathbb{R}^d &\to \mathbb{R}^{N-1}\\
x&\mapsto \Psi_t(x) = \Big(\lambda_1^t\psi_1(x),
\dots,\lambda_{N-1}^t\psi_{N-1}(x)\Big)^\intercal,
\end{align*}
which can be interpreted as a feature map embedding the data points in $\mathbb{R}^{N-1}$. The fundamental property relating the diffusion maps to the diffusion distance is that the diffusion distance is the Euclidean distance in the embedding space $\mathbb{R}^{N-1}$, i.e.,
 \[
D^2_t(x,z) =  \|\Psi_t(x)-\Psi_t(z)\|_2^2,
\]
which can be intuitively understood as an ``isometry'' property, as it was mentioned in the introduction of this paper.
Notice that we excluded here the term $\psi_0(x) =1$ which is constant and, therefore, does not influence the diffusion distance. Indeed, the diffusion map is an embedding of each data point into $\mathbb{R}^{N-1}$, so that, to put it in a nutshell, two data points are close to each other if the time to diffuse from one to the other is small.

The natural approach consists in the approximation $\lambda_\ell^t\approx 0$ for all $\ell> q(t)$, so that the diffusion maps become
\begin{equation}
x\mapsto \Big(
\lambda_1^t\psi_1(x),\dots,\lambda_{q(t)}^t\psi_{q(t)}(x)\Big)^\intercal,\label{eq:TruncatedDiffusionMap}
\end{equation}
with $q(t)\ll N$, where we have an orthogonality relation of the components with respect to a non-canonical inner product, i.e., $\langle\psi_\ell,\psi_{\ell'}\rangle_{\phi^{-1}} = 0$, for all $\ell\neq \ell'$. The spectral truncation of the transition matrix~(\ref{eq:Transition Proba})
\begin{equation}
p^\star_{t}(x,y) =  \sum_{\ell =0}^{m-1}\lambda_{\ell}^t \psi_{\ell}(x)\phi_\ell(y),\label{eq:OptimalTruncation}
\end{equation}
is known to be an optimal approximation in $m$ dimensions in the least squares sense (see, e.g.,~\cite{Nadler}).
\begin{mythm}[Optimality of spectral truncation~\cite{Nadler}]\label{Thm1}
The truncated transition matrix~(\ref{eq:OptimalTruncation}) is the solution to the least squares minimization problem
\[
\min_{\hat{p}}\mathbb{E}_{x}\Big\{\|p_{t}(x,\cdot)-\hat{p}_{t}(x,\cdot) \|^2_{\phi_0^{-1}}\Big\},
\]
over $m$-dimensional approximations of the form $\hat{p}_{t}(x,y) =  \sum_{\ell =0}^{m-1}a_{\ell}(x)w_\ell(y)$, where the expectation is taken with respect to the stationary probability density $\phi_0(x)$.
\end{mythm}
An important remark about spectral truncation is that it is in fact solving a minimization of a weighted sum of squares, where the weight is given by $\phi_0(x)$, i.e., the weight is proportional to the volume of the data points. Indeed, the spectral truncation will generally speaking make larger errors on data point with a small volume (or small stationary probability) which can possibly be outliers. Hence, the spectral truncation will be more accurate in dense regions where the diffusion process has a large probability at stationarity.

Let us notice that the inner product between the position vectors of $x$ and $y$  in the embedding space at time $t$ is related to the transition probability at time $t$ as follows
\begin{equation}
\Psi^\intercal_t(x)\Psi_t(y) = \sum_{\ell =1}^{N-1}\lambda_{\ell}^{2t} \psi_{\ell}(x)\psi_\ell(y)= p_{2t}(x,y)\phi^{-1}_0(y)-1,\label{eq:DiffKernel}
\end{equation}
where the term $\ell =0$ (with $\lambda_0=1$) is absent.
By analogy with kernel methods, we define the ``diffusion kernel'' matrix $
K_t(x,y) = \Psi^\intercal_t(x)\Psi_t(y)$,
which is a positive semi-definite matrix of $x,y\in \Omega$. Hence, the diffusion map is a feature map related to the diffusion kernel. In the context of kernel methods, the squared diffusion distance would be called  squared ``kernel distance'' between two data points $x$ and $y$, i.e.,
 \[
D^2_t(x,y) =  K_t(x,x)+K_t(y,y)-2K_t(x,y).
\]
The diagonal element $K_t(x,x)$ is related to the probability to come back at position $x$ at time $t$ while starting at $x$ at time $0$. In particular for $t=1/2$, the length of the position vector $x$ in the embedding space   \[\|\Psi_{1/2}(x)\|_2 = \frac{1}{d_\Omega(x)}\Big(k(x,x)\vol( \Omega)-d^2_\Omega(x)\Big)^{1/2}\] is, in approximation, inversely proportional to the volume of $x$. Intuitively, this means that a point with a large probability of presence at stationarity (a ``central'' data point) will be embedded close to $0$ in the embedding space.

A non-trivial result relating the squared degree of each data point to the volume of the whole dataset is easily obtained as an indirect consequence of the positive definiteness of the kernel $k$.
\begin{myprop}\label{Prop:2}
Let $x\in \Omega$ and let $k$ be a Mercer kernel such that $k(x,y)>0$, for all $x,y\in X$. Then, it holds that $d^2_\Omega(x)\leq k(x,x)\vol( \Omega)$,
or more explicitly, we have the inequality
 \[
\sum_{y,y'\in \Omega}k(x,y)k(y',x)\leq k(x,x)\sum_{y,y'\in \Omega}k(y,y').
\]

\end{myprop}
\begin{proof}
From the identity~(\ref{eq:DiffKernel}), we find
$d^2_\Omega(x)\|\Psi_{1/2}(x)\|^2_2 = k(x,x)\vol( \Omega)/d^2_\Omega(x)-1\geq 0$,
which merely yields $d^2_\Omega(x)\leq k(x,x)\vol( \Omega)$.
\end{proof}
Actually, Proposition~\ref{Prop:2} will be useful in order to define a proper extension of the kernel $K_{1/2}$ to new data points in $X\setminus \Omega$ as it is explained in Section~\ref{sec:oos}. Let us now recall some aspects of the duality theory of SDP's.

\section{Useful properties of the solution of the semi-definite program\label{sec:SDP}}
In this section, several results are consequences of the papers~\cite{Alizadeh1997,Burer2003}.
In order to have a different interpretation, as detailed in~\ref{App:Dual},  we also formulate  the Lagrange dual of~(\ref{eq:OptimalityProperty}), that we recall here
\[\max_{\rho \succeq 0}\Tr (\rho K) ,\ {\rm subject\ to}\ \rho(x,x) =K(x,x)\ {\rm for\ all}\ x\in \Omega. \]
This Lagrange dual reads
\begin{equation}
\min_{L\succeq 0} \Tr \big(\ddiag(K)L\big)+\Tr \big(\ddiag(K)K\big), \ {\rm s.t.}\  L+ K \text{ is diagonal}.
\label{eq:DualFixedOnesDiagonalApproximation}
\end{equation}
By the diagonality constraint, we mean $[L+K]_{ij} = 0$ for all $i\neq j$ in $\{1,\dots, N\}$.
Indeed, the dual optimization problem aims to obtain the \emph{p.s.d.} matrix of least trace which only differs from $-K$ by its diagonal elements.  Following~\cite{Alizadeh1997}, the primal problem~(\ref{eq:OptimalityProperty}) is strictly feasible since $\rho = \ddiag(K)$ is feasible and strictly positive definite because the diagonal elements are all strictly positive reals. Notice that the dual problem is also strictly feasible~\cite{Bandeira2016}.
The connection between the primal and dual solutions is given by the following result.
\begin{myprop}[Dual certificate]\label{Prop:DualCertif}
A \emph{p.s.d.} symmetric matrix $\rho^\star\in \mathbb{R}^{N\times N}$ satisfying
$\rho^\star(x,x) = K(x,x)$,  for all $x\in \Omega$,
is a global optimum of (\ref{eq:OptimalityProperty}) if and only if there exists a symmetric $L^\star\in \mathbb{R}^{N\times N}$ such that:
\begin{enumerate}
\item[(i)] $L^\star\rho^\star = 0$,
\item[(ii)] $L^\star\succeq 0$,
\item[(iii)] $L^\star+K$ is diagonal.
\end{enumerate}
\end{myprop}
In~\cite{Bandeira2016,BoumalGeNPowMeth}, $L^\star$ is called dual certificate. The notation $L$ is used in order to emphasize that this matrix shares similarities with a Laplacian matrix as noticed in~\cite{Bandeira2016}, with the difference that $K$ is not a weight matrix. As explained generally in~\cite{Alizadeh1997}, because of the complementary slackness condition, $L^\star$ and $\rho^\star$ commute and therefore they share the same basis of eigenvectors.
As a simple consequence of $L^\star\rho^\star = 0$, the eigenvectors of non-zero eigenvalues of $\rho$  are a basis of ${\rm ker} L^\star$. Hence, these eigenvectors play an analogous role to the indicator functions of the connected components of a graph in the context of spectral graph theory.

To summarize, we have a unique solution for both the primal and dual problems. There is no duality gap, i.e., the optimal objective values are equal: $ \Tr \big({\ddiag(K)}L^\star\big)+\Tr \big({\ddiag(K)}K\big)= \Tr(K\rho^\star)$.
Following the same argument as in~\cite{BoumalGeNPowMeth}, we can show that the diagonal matrix $D = L^\star + K$ is,
\[
D(x,x) = K(x,x)^{-1}(K\rho^\star)(x,x),  \ {\rm for\ all}\ x\in \Omega.
\]
In view of applications, our objective is to use these properties in the context of diffusion geometry.

The previous results gave a necessary and sufficient condition for a \emph{p.s.d.} matrix to be optimal. However, the solution of the SDP~(\ref{eq:OptimalityProperty}) is not always interesting as announced in Remark~\ref{Rem:Positivity}. Indeed, an elementary example concerns the specific choice of a kernel taking only non-negative values, i.e., $
K(x,y)\geq 0$,  for all $x,y\in\Omega$. Indeed, this choice always yields a trivial rank one optimal solution $\rho^\star$.

\begin{mylemma}[Trivial solution]\label{Lem:Trivial} Let $\rho_K$ be the rank one symmetric and \emph{p.s.d.} matrix given by $\rho_K(x,y) = \sqrt{K(x,x)}\sqrt{K(y,y)}$ for all $x,y\in\Omega$.
Then, we have
\begin{enumerate}
\item[(i)] $L(\rho_K)\rho_K = 0$,
\item[(ii)] and if the kernel $K$ further satisfies $K(x,y)\geq 0$ for all $x,y\in\Omega$, then $L(\rho_K) \succeq 0$,
\end{enumerate}
where $L(\rho) = \ddiag(K)^{-1}\ddiag(K\rho)-K$ is the matrix defined in~(\ref{eq:DualCertificate}).
\end{mylemma}
\begin{proof}Let $x,y\in\Omega$. We firstly have
\[
L(\rho_K)(x,y) = \delta_{x,y}\sum_{z\in\Omega}K(x,z)\sqrt{\frac{K(z,z)}{K(x,x)}}-K(x,y),
\]
where $\delta_{x,y}$  is the Kronecker symbol.
In order to prove the first property, we simply have that
\[
\sum_{y\in\Omega}L(\rho_K)(x,y)\sqrt{K(y,y)} =0,
\]
from the definition of $L(\rho_K)$. In order to prove the second property, let $v:\Omega\to \mathbb{R}$. Then, we simply verify that
\[
\sum_{x,y\in\Omega}v(x)L(\rho_K)(x,y)v(y) = \frac{1}{2}\sum_{x,y\in\Omega}K(x,y)\Big(\sqrt{\frac{K(y,y)}{K(x,x)}}v(x)-\sqrt{\frac{K(x,x)}{K(y,y)}}v(y)\Big)^2\geq 0,
\]
if $K(x,y)\geq 0$ for all $x,y\in\Omega$.
\end{proof}
Hence, the choice of a kernel yielding a non trivial solution is restricted to kernels taking positive as well as negative values.
Lemma~\ref{Lem:Trivial} provides an example where the optimal solution has rank $1$. Let us now comment about the rank of the optimal solution $\rho^\star$ in another circumstance. In view of applications, we make the specific choice
\[
K(x,y) =\frac{k(x_,y)}{\sqrt{d_\Omega(x)}\sqrt{d_\Omega(y)}}-\sqrt{\frac{d_\Omega(x)}{\vol(\Omega)}}\sqrt{\frac{d_\Omega(y)}{\vol(\Omega)}},
\]
which is positive semi-definite with a maximal eigenvalue strictly smaller than one. Under these assumptions, Proposition~\ref{Prop:Rank} states that finding $\rho^\star$ is in fact equivalent to solving a convex relaxation of a rank minimization problem in terms of the nuclear norm $\|\cdot\|_{\star}$. Indeed, since the nuclear norm is a typical convexification of the rank, Proposition~\ref{Prop:Rank} illustrates why in practice the optimal solution of~(\ref{eq:OptimalityProperty}) is of low rank.
\begin{myprop}[Equivalence with a nuclear norm minimization]\label{Prop:Rank}
Let $K\in \mathbb{R}^{N\times N}$ be a \emph{p.s.d.} matrix with a maximal eigenvalue strictly smaller than $1$ and let $\Sigma \in \mathbb{R}^{N\times N}$ be the invertible matrix with orthogonal rows such that $\mathbb{I}-K= \Sigma\Sigma^\intercal$. Then, the optimal solution $X^\star$ of
\[
\min_{X\succeq 0}\|X\|_{\star} ,\ {\rm subject\ to}\ \ddiag\Big( (\Sigma^{-1})^\intercal X\Sigma^{-1}\Big) = \ddiag(K),
\]
has the same rank as the optimal solution of~(\ref{eq:OptimalityProperty}) $\rho^\star$ and is given by $X^\star = \Sigma^\intercal \rho^\star\Sigma$.
\end{myprop}
\begin{proof}
Notice first that maximizing $\Tr(\rho K)$ over $\rho\succeq 0$ such that $\ddiag(\rho) = \ddiag(K)$ is equivalent to minimizing $\Tr\big(\rho(\mathbb{I}- K)\big)$ since $\Tr(\rho)$ is fixed by the constraints. Then, we remark that $(\mathbb{I}- K)\succ 0$ because the maximal eigenvalue of $K$ is strictly smaller than $1$ by assumption. Hence, a diagonalization procedure yields $\mathbb{I}- K = U D U^\intercal$ where the diagonal matrix $D$ is strictly positive definite. Therefore, we can also write $\mathbb{I}- K = \Sigma\Sigma^\intercal$,  where $\Sigma\in \mathbb{R}^{N\times N}$ is invertible. As a consequence, the objective of the minimization problem
\[
\min_{\rho\succeq 0} \Tr\big(\rho(\mathbb{I}- K)\big),\ {\rm subject\ to}\ \ddiag(\rho)= \ddiag(K),
\]
can be written $\Tr\big(\rho(\mathbb{I}- K)\big) = \Tr\big(\Sigma^\intercal\rho\Sigma\big)$, where $\Sigma^\intercal\rho\Sigma\succeq 0$.  Then, by performing the change of variables $X = \Sigma^\intercal\rho\Sigma$, we can rephrase the minimization problem as follows
\[
\min_{X\succeq 0} \Tr\big(X),\ {\rm subject\ to}\ \ddiag\Big((\Sigma^{-1})^\intercal X\Sigma^{-1}\Big)= \ddiag(K).
\]
Finally, we notice that $\|X\|_\star = \Tr(\sqrt{XX^\intercal}) = \Tr(X)$ since $X\succeq 0$.
\end{proof}

Importantly, solving directly the nuclear norm minimization of Proposition~\ref{Prop:Rank} is not feasible in practice for large scale matrices since it requires a full diagonalization of $K$.
\begin{myremark}[Analogy with quantum mechanics]
The objective function in~(\ref{eq:OptimalityProperty}) can be bounded as follows
$\Tr(\rho K) \leq \sqrt{\Tr(\rho^2)} \sqrt{\Tr(K^2)}$.
Intuitively, the maximization~(\ref{eq:OptimalityProperty}) can yield solutions with a large value of the Frobenius norm $\sqrt{\Tr(\rho^2)}$. Since the constraints fix the diagonal of $\rho$, the trace is fixed to $\Tr(\rho) = \Tr(K)$. Then, the normalized matrix
\[
\hat{\rho} = \rho/\Tr(\rho) = \rho /\Tr(K),
\] can be equivalently considered. In quantum mechanics, $\hat{\rho}$ is named a density matrix and $\Tr(\hat{\rho}^2)$ is called the ``purity'' of the density matrix $\hat{\rho}$. Typically, $\rho$ is then a complex Hermitian p.s.d. matrix. A pure state corresponds to the rank $1$ case $\Tr(\hat{\rho})=\Tr(\hat{\rho}^2)$. In this context, fixing the diagonal of the density matrix amounts to fix the probability distribution of the mixed state represented by $\hat{\rho}$. Finally, the objective function $\Tr (\hat{\rho} K)$ is the expectation value of the observable $K$ for the quantum system described by $\hat{\rho}$.
\end{myremark}

\section{Out-of-sample extension of the normalized diffusion Kernel \label{sec:oos}}
Let us discuss now the out-of-sample extension of data points for the diffusion kernel $K$ defined in~(\ref{eq:DiffKernel}) and associated to the discrete set $\Omega=\{x_1,\dots,x_N\}\subset X$.
Let us recall that the diffusion kernel $K$ defined on $\Omega\times \Omega$ is given by  \[
K(x_i,x_j) = k_\Omega(x_i,x_j)- e_0(x_i)e_0(x_j),\quad \text{ for all } x_i,x_j\in \Omega,
\]
where the normalized kernel is
\[
k_\Omega(x_i,x_j) = \frac{k(x_i,x_j)}{\sqrt{d_\Omega(x_i)}\sqrt{d_\Omega(x_j)}}
\]
and with $e_0(x)$ is the normalized eigenvector with eigenvalue $1$ of $k_\Omega(x,y)$ defined on $\Omega\times \Omega$, given explicitly by $
e_0(x_i) = \sqrt{d_\Omega(x_i)/\vol(\Omega)}$. The diagonal elements $(x_i,x_i)$ of the kernel
$K(x_i,x_i) = 1/d_\Omega(x_i)-d_\Omega(x_i)/\vol(\Omega)$,
is the difference between the probability that the lazy random walker does not jump at time $t+1$ (i.e., remains at $x_i$) with the probability to find the random walker at $x_i$ after an asymptotically long time. Indeed, the first term only depends locally of the neighborhood of $x_i$ since the Gaussian kernel decays off very fast, while the second term includes the influence of all the points in the whole data set.  As a consequence of Proposition~\ref{Prop:2}, this difference is always a positive number.

Let us now consider the extension of the kernel $K$ to points in $X$. Firstly, it is straightforward to introduce the extension of the degree function as follows \[\bar{d}_\Omega(x)=\sum_{z\in \Omega}k(x,z) \quad \text{ for all } x\in X,\]
by using the kernel $k$ naturally defined on $X\times X$. Similarly, the extension of the normalized kernel is straightforward and was already proposed in~\cite{Bengio}. Indeed, we have $\overline{k_\Omega}:X\times X\to\mathbb{R}$ such that
\[
\overline{k_\Omega}(x,y) = \frac{k(x,y)}{\sqrt{\bar{d}_\Omega(x)}\sqrt{\bar{d}_\Omega(y)}} ,
\]
for all $x,y\in X$. Let us emphasize that the extension of $k_\Omega$ clearly yields a \emph{p.s.d.} kernel function and in that sense, it is particularly natural. However, as it was analysed in the previous section, the object of interest is here the diffusion kernel $K_{1/2}$ whose extension is less trivial. In order to define it, we firstly deal with the extension of the top eigenvector of $k_\Omega$. The extension of the eigenvector $\overline{e_0}:X\to \mathbb{R}$ is then given by the Nystr\"om formula~\cite{WilliamSeeger}
\[
\overline{e_0}(x) = \sum_{z\in \Omega}\overline{k_\Omega}(x,z)e_0(z)= \sqrt{\frac{\bar{d}_\Omega(x)}{\vol(\Omega)}}\quad \text{ for all } x\in X.
\]
Notice that $\overline{e_0}(x)$ is in general not an eigenvector of an integral operator associated to $\overline{k_\Omega}:X\times X\to\mathbb{R}$.
However, the following extension of the diffusion kernel can be defined
$
\overline{K}(x,y) = \overline{k_\Omega}(x,y)- \overline{e_0}(x)\overline{e_0}(y)
$
for all $x,y\in X$, and in particular the diagonal of the extended kernel reads
\begin{equation}\label{eq:diagonal_extension}
\overline{K}(x,x) = \frac{k(x,x)}{\bar{d}_\Omega(x)}-\frac{\bar{d}_\Omega(x)}{\sum_{z\in \Omega}\bar{d}_\Omega(z)}, \quad x\in X,
\end{equation}
and where $k(x,x) =1$ in the case of the Gaussian kernel.
Since the diffusion kernel is intrinsically normalized with respect to the data set $\Omega$ so that it is associated to a diffusion process on $\Omega$, its out-of-sample extension intuitively makes sense if the volume of the out-of-sample point is small enough with respect to the total volume of the data set $\Omega$. Indeed, the extension formula~(\ref{eq:Out-of-sample2}) requires that $\bar{K}(x,x)\geq 0$ for all $x\in X$ so that its square root can remain real. This requirement is \emph{a priori} not trivially satisfied.
A bit surprisingly, Proposition~\ref{Prop:3} states that the extension~(\ref{eq:diagonal_extension}) of the diagonal of the kernel is always well-defined, that is, the volume of any new data point is always small enough with respect to the volume of the dataset.
\begin{myprop}\label{Prop:3}
Let $\bar{x}\in X$. Then, we have
$\overline{K}(\bar{x},\bar{x})\geq 0$. Equivalently, it holds that $\big(\bar{d}_\Omega(\bar{x})\big)^2\leq k(\bar{x},\bar{x}) \vol( \Omega)$.
\end{myprop}
\begin{proof}
Let us consider the same diffusion process as before but defined on $\bar{\Omega} =\Omega\cup \{\bar{x}\}$.  We firstly notice that the connection between the squared degrees \[d^2_{\bar{\Omega}}(\bar{x}) = \Big(k(\bar{x},\bar{x})+\sum_{y\in \Omega}k(\bar{x},y)\Big)^2 =k^2(\bar{x},\bar{x})+2k(\bar{x},\bar{x})\bar{d}_{\Omega}(\bar{x})+ \big(\bar{d}_{\Omega}(\bar{x})\big)^2,\]
and the relation between the total volumes,
\[\vol( \bar{\Omega}) = d_{\bar{\Omega}}(\bar{x})+\sum_{y\in \Omega}d_{\bar{\Omega}}(y)  = k(\bar{x},\bar{x})+\sum_{y\in \Omega} k(y,\bar{x})+\sum_{y\in \Omega}\Big(\sum_{x\in\Omega} k(y,x) + k(y,\bar{x}) \Big)  = k(\bar{x},\bar{x})+2\bar{d}_{\Omega}(\bar{x})+\vol(\Omega).\]
Then, as a consequence of Proposition~\ref{Prop:2}, we have $d^2_{\bar{\Omega}}(\bar{x})\leq k(\bar{x},\bar{x})\vol( \bar{\Omega})$, which implies after a simplification the desired result: $
 \big(\bar{d}_{\Omega}(\bar{x})\big)^2\leq k(\bar{x},\bar{x})\vol(\Omega)$.
\end{proof}
Finally, we mention that although $\overline{k_\Omega}(x,y)$ is a Mercer kernel, $\overline{K}(x,y)$ does not have this property. This does not cause any difficulty for the out-of-sample formula~(\ref{eq:Out-of-sample2}) since only the diagonal of $\overline{K}(x,y)$ has to remain positive.
\begin{myremark}
Notice that, in general, the extended diffusion kernel
\[
\overline{K}(x,y) = \frac{k(x,y)}{\sqrt{\bar{d}_\Omega(x)\bar{d}_\Omega(y)}}-\frac{\sqrt{\bar{d}_\Omega(x)\bar{d}_\Omega(y)}}{\sum_{z\in \Omega}\bar{d}_\Omega(z)}, \quad x, y\in X,
\]
is indefinite although its diagonal is positive. This is seen by choosing to evaluate this kernel on $\Omega\cup\{\bar{x}\}$. Indeed, we have
\[
\sum_{x,y\in \Omega\cup\{\bar{x}\}}\sqrt{\bar{d}_\Omega(x)}\overline{K}(x,y)\sqrt{\bar{d}_\Omega(y)} = \sum_{y\in \Omega\cup\{\bar{x}\}}\bar{d}_\Omega(x)\Big(1-\frac{\sum_{x\in \Omega\cup \{\bar{x}\}}\bar{d}_\Omega(x)}{\sum_{z\in \Omega}\bar{d}_\Omega(z)}\Big)\leq 0,
\]
since $\sum_{x\in \Omega\cup \{\bar{x}\}}\bar{d}_\Omega(x)\geq \sum_{z\in \Omega}\bar{d}_\Omega(z)$.
\end{myremark}
We now address the numerical solution of the optimization problem~(\ref{eq:OptimalityProperty}).
\section{Simulations\label{sec:NumStrat}}

The outline of our method is as follows. The kernel $K$ as it is given in~(\ref{eq:KDiff}) is used in order to define the semi-definite program~(\ref{eq:OptimalityProperty}).
Then, the embedding coordinates are calculated by computing the eigenvectors of $\rho^\star$. Finally, the out-of-sample formula can be used thanks to the result of the previous section.

\subsection{Numerical method}
In order to solve numerically the problem~(\ref{eq:OptimalityProperty}), a reformulation is firstly proposed. Namely, we define the positive semi-definite matrix \[\tilde{\rho} \triangleq \ddiag(K)^{-1/2}\rho\ddiag(K)^{-1/2}\] which is supposed  to be decomposed as follows $\tilde{\rho} = HH^\intercal$, with $H\in \mathbb{R}^{N\times r_0}$. We further define a coupling matrix $J = \ddiag(K)^{1/2}K\ddiag(K)^{1/2}$, as well as a quadratic objective function $E(H) = \Tr (H^\intercal JH)$. Under the assumption that $\tilde{\rho} = HH^\intercal$, we can equivalently solve
\begin{equation}
\max_{H\in \mathbb{R}^{N\times r}}E(H)\ {\rm s.t.} \ \sum_{\ell=1}^{r}[H]_{i\ell}^2 = 1, \text{ for all } i\in\{1,\dots,N\}\label{eq:factorizedProblem}
\end{equation}
yielding a solution $H^\star$ related to the candidate optimal solution of~(\ref{eq:OptimalityProperty}) by the relation \[\rho^\star = \ddiag(K)^{1/2}H^\star {H^\star}^\intercal \ddiag(K)^{1/2}.\]
For convenience, let us write $H_\Xi = \ddiag(K)^{1/2}H^\star$.
Indeed, the objective function $E(H)$ is quadratic in $H$ while the feasible space is a product of $N$ unit spheres.
Inspired by the projected gradient method, a natural projection operator from $\mathbb{R}^{N\times r_0}$ to the feasible $\mathcal{M} = (\mathcal{S}^{r_0-1})^N$ is simply obtained by projecting each factor of the Cartesian product on the unit sphere $ \mathbb{S}^{r_0-1}$, i.e., let $i \in\{1,\dots, N\}$, then
\begin{equation}
[\mathcal{P}(H)]_{i\ell} = [H]_{i\ell}/\Big(\sum_{\ell'=1}^{r_0}[H]_{j\ell'}^2\Big)^{1/2},\label{eq:P}
\end{equation}
normalizes the rows of the matrix $H\in \mathbb{R}^{N\times r_0}$. Hence, the method in order to maximize $E(H)$ summarized in Algorithm~\ref{Alg1} consists of a succession of matrix multiplications by $J$ and projection steps $\mathcal{P}$.
\begin{algorithm}[h]
\caption{Projected power method~\cite{BoumalGeNPowMeth} \label{Alg1}}
\begin{algorithmic}[1]
\Require Symmetric positive definite matrix $J\in \mathbb{R}^{N\times N}$;
and an initial $H_0\in \mathbb{R}^{N\times r_0}$ such that $\mathcal{P}(H_0) = H_0$.
\State {\bf for} $n =1, 2, \dots$ {\bf do}
\State $H_n = \mathcal{P}(JH_{n-1})$.
\State  {\bf end for}
\end{algorithmic}
\end{algorithm}
The sequence of iterates of Algorithm~\ref{Alg1} have increasing objective values, as explained in~\cite{BoumalGeNPowMeth}. Once a solution $H^\star$ is found by using Algorithm~\ref{Alg1}, the embedding coordinates are found by computing the singular value decomposition of $H_\Xi$.

To summarize, the numerical algorithm used to solved a rank constrained version of the semi-definite program~(\ref{eq:OptimalityProperty})  is the following:
\begin{enumerate}
\item A matrix $M_0\in \mathbb{R}^{N\times r_0}$ is generated with independent entries in $[-1,1]$ chosen uniformly at random. Then, the initial point for Algorithm~\ref{Alg1} is obtained as $H_0 = \mathcal{P}(M_0)$.
\item The projected power method is applied to a  problem yielding a result $H^\star \in \mathbb{R}^{N\times r_0}$ after convergence.
\item Then, the optimality of the candidate solution $\rho^\star = H_\Xi H_\Xi^\intercal$ with $H_\Xi = \ddiag(K)^{1/2}H^\star$ is verified by calculating $L(\rho^\star) H_\Xi \in \mathbb{R}^{N\times r_0}$ and the $6$ least eigenvalues of $L(\rho^\star)$.
\item Finally, a singular value decomposition of $H_\Xi$ is performed in order to obtain the embedding coordinates.
\end{enumerate}
Incidentally, in the context of graph embedding, the maximization problem~(\ref{eq:factorizedProblem}) has similarities  with a relaxation of modularity maximization studied in the master thesis~\cite{MasterThesis}.

\subsection{Illustration of the optimal kernel\label{sec:Toy}}
We choose here to illustrate the problem~(\ref{eq:FormalDefSDP}) on a simple toy model, while a reader interested in applications of dimensionality reduction is referred to the next subsection. We choose here $X = [-1,1]$ with the Lebesgue measure. Then, the initial continuous positive semi-definite kernel is defined as follows:
\[
K(x,y) = \frac{e^{-(x-y)^2/\sigma^2}}{\sqrt{d(x)d(y)}}-\frac{\sqrt{d(x)d(y)}}{\iint_{-1}^1e^{-(x-y)^2/\sigma^2}\mathrm{d} x\mathrm{d} y},\text{ with } d(x) = \int_{-1}^1 e^{-(x-y)^2/\sigma^2}\mathrm{d} y
\]
 and $\sigma >0$. Noticeably, in accordance with Remark~\ref{Rem:Positivity}, $K(x,y)$ is chosen to take both positive and negative value. Then, the square $[-1,1]\times[-1,1]$ is discretized in a square grid of $2001\times 2001$ points. The result of the discrete optimization problem is displayed in Figure~\ref{Fig:Toy1} for $\sigma = 0.1$ and in Figure~\ref{Fig:Toy2} for $\sigma = 0.7$. The eigenvectors of the rank two solution $\rho^\star$ are compared to the two leading eigenvectors of $K$ in Figure~\ref{Fig:Toy_eigs1} and  Figure~\ref{Fig:Toy_eigs2}, for $\sigma = 0.1$ and $\sigma = 0.7$, respectively. It is interesting to notice in that case the qualitative similarity of the eigenfunctions. In Figure~\ref{Fig:EmbedInter01} and Figure~\ref{Fig:EmbedInter07}, the discretized interval is embedded in the plane by using these eigenfunctions.

\begin{figure}[h]
\centering
\includegraphics[trim={2cm 8.5cm 2cm 8.5cm},clip,scale = 0.45]{Toy_K_3D_red.pdf}
\includegraphics[trim={2cm 8.5cm 2cm 8.5cm},clip,scale = 0.45]{Toy_rho_3D_red.pdf}
\caption{For $\sigma = 0.1$, on the left, $K(x,y)$ and on the right $\rho^\star(x,y)$ are plotted on the square $[-1,1]\times[-1,1]$.\label{Fig:Toy1}}
\end{figure}
Let us emphasize that the diagonal of the kernels $K$ and $\rho^\star$ is preserved by the numerical approach, which uses the procedure described in the previous section.
\begin{figure}[h]
\centering
\setlength\figureheight{0.15\textwidth}
\setlength\figurewidth{0.3\textwidth}
\begin{minipage}{0.49\textwidth}
\input{Eigs_Toy_K.tikz}
\end{minipage}
\begin{minipage}{0.49\textwidth}
\input{Eigs_Toy_rho.tikz}
\end{minipage}
\caption{For $\sigma = 0.1$, on the left, the two leading eigenvectors of the discretized  $L_K$ and on the right, the two eigenvectors of non-zero eigenvalues of the discretized $L_{\rho^\star}$ (i.e., $\rank(\rho^\star) = 2$). The blue curve represents the eigenvector associated with the largest eigenvalue ($e_0$ or $\chi_1$) and is an odd function, whereas the red curve is an even function ($e_1$ or $\chi_2$).\label{Fig:Toy_eigs1}}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[trim={2cm 8.5cm 2cm 8.5cm},clip,scale = 0.45]{SegmentDM01.pdf}
\includegraphics[trim={2cm 8.5cm 2cm 8.5cm},clip,scale = 0.45]{SegmentSDPEmbedding01.pdf}
\caption{Embedding of the discretized interval $[-1,1]$ for $\sigma = 0.1$ thanks to Diffusion Maps (lhs) and the SDP embedding (rhs).\label{Fig:EmbedInter01}}
\end{figure}
In Figure~\ref{Fig:Toy1}, the numerical solution obtained for $\rho^\star (x,y)$ is a continuous function on $[-1,1]\times[-1,1]$. In contrast with the case $\sigma =0.1$, we observe in Figure~\ref{Fig:Toy2} that the numerical solution obtained for $\rho^\star (x,y)$ exhibits a step-like structure which is close to be a discontinuity, for the bandwidth $\sigma = 0.7$. If we increase the bandwidth to $\sigma = 1$, we obtain numerically the optimal solution \[
\rho^\star(x,y) = \sign(e_0(x))\sqrt{K(x,x)}\sign(e_0(y))\sqrt{K(y,y)},
\]
which is of rank one and where $\sign(e_0(x)) = \sign(x)$ for all $x\in[-1,1]$ (we recall that $e_0$ is the eigenfunction of $L_K$ with the largest eigenvalue). Its form can be guessed from the shape of the eigenvectors in Figure~\ref{Fig:Toy_eigs2}. Noticeably, this kernel saturates the first upper bound in~(\ref{eq:UniBounded}), i.e.,
$|\rho(x,y)| \leq\sqrt{K(x,x)}\sqrt{K(y,y)}$,
which is a consequence of the positive definite chararacter of $\rho$.
\begin{figure}[H]
\centering
\includegraphics[trim={2cm 8.5cm 2cm 8.5cm},clip,scale = 0.45]{Toy_K_3D_07_red.pdf}
\includegraphics[trim={2cm 8.5cm 2cm 8.5cm},clip,scale = 0.45]{Toy_rho_3D_07_red.pdf}
\caption{For $\sigma = 0.7$, on the left, $K(x,y)$ and on the right $\rho^\star(x,y)$ are plotted on the square $[-1,1]\times[-1,1]$.\label{Fig:Toy2}}
\end{figure}
\begin{figure}[H]
\centering
\setlength\figureheight{0.15\textwidth}
\setlength\figurewidth{0.3\textwidth}
\begin{minipage}{0.49\textwidth}
\input{Eigs_Toy_K_07.tikz}
\end{minipage}
\hfill
\begin{minipage}{0.49\textwidth}
\input{Eigs_Toy_rho_07.tikz}
\end{minipage}
\caption{For $\sigma = 0.7$, on the left, the two leading eigenvectors of the discretized  $L_K$ and on the right, the two eigenvectors of non-zero eigenvalues of the discretized $L_{\rho^\star}$. The blue curve represents the eigenvector associated with the largest eigenvalue ($e_0$ or $\chi_1$) and is an odd function, whereas the red curve is an even function ($e_1$ or $\chi_2$).\label{Fig:Toy_eigs2}}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[trim={2cm 8.5cm 2cm 8.5cm},clip,scale = 0.45]{SegmentDM07.pdf}
\includegraphics[trim={2cm 8.5cm 2cm 8.5cm},clip,scale = 0.45]{SegmentSDPEmbedding07.pdf}
\caption{Embedding of the discretized interval $[-1,1]$ for $\sigma = 0.7$ thanks to Diffusion Maps (lhs) and the SDP embedding (rhs).\label{Fig:EmbedInter07}}
\end{figure}

\subsection{Illustrative examples of dimensionality reduction\label{sec:Data}}
By using the numerical strategy described hereabove, we illustrate here the embedding method on several real and artificial data. Let us emphasize that the only parameter to be selected is the bandwidth parameter $\sigma$ of the Gaussian kernel $k(x,y) = \exp(-\|x-y\|_2^2/\sigma^2)$ used to defined the diffusion kernel~(\ref{eq:DiffKernel}).

\begin{figure}[H]
\setlength\figureheight{0.2\textwidth}
\setlength\figurewidth{0.3\textwidth}
\begin{minipage}{0.49\textwidth}
\input{originaldataoutliers.tikz}
\end{minipage}\hfill
\begin{minipage}{0.49\textwidth}
\input{SpectrumDMdataoutliers.tikz}
\end{minipage}
\begin{minipage}{0.49\textwidth}
\input{SDPdataoutliers.tikz}
\end{minipage}\hfill
\begin{minipage}{0.49\textwidth}
\input{DMdataoutliers.tikz}
\end{minipage}
\begin{minipage}{0.49\textwidth}
\input{SDPdataoutlierssigma5.tikz}
\end{minipage}\hfill
\begin{minipage}{0.49\textwidth}
\input{DMdataoutlierssigma5.tikz}
\end{minipage}
\caption{Illustration of the rigidity of the SDP embedding for several values of $\sigma$ with respect to Diffusion Maps. The dataset contains three clusters with 100 points and $8$ outliers (blue stars). The original data is on top left. The spectrum of the kernel~(\ref{eq:DiffKernel}) on the top right estimates the quality of the spectral truncation of the diffusion maps. The color is the probability of presence at stationarity $\phi_0(x)$. In all cases, we have $\rank(\rho^\star) = 2$ for the SDP embedding.\label{Fig:OutliersBandwidth}}
\end{figure}
The datasets  used in this section include data points both with homogeneous and heterogeneous attributes. Let us shortly describe them:
\begin{itemize}
\item Iris: 150 points in dimension $4$ (\url{https://archive.ics.uci.edu/ml/datasets/iris}).
\item Wine: 178 points in dimension 13 (\url{https://archive.ics.uci.edu/ml/datasets/wine}).
\item MNIST: pictures of hand written digits, corresponding to dimension 784 (\url{http://yann.lecun.com/exdb/mnist/}).
\end{itemize}
\begin{figure}[H]
\centering
\begin{tikzpicture}
\draw (-2,1)node[right]{\bf Input};
\draw (-2,0.5) node[right]{$\Omega =\{x_1,\dots,x_N\}\subset X$};
\draw (-2,0) node[right]{$K:X\times X\to \mathbb{R}$ cfr.~(\ref{eq:KDiff})};
\draw (-2,-0.5) node[right]{ and $K\in\mathbb{R}^{N\times N}$ p.s.d.};
\draw (-2,-1) node[right]{kernel matrix.};
\draw[thick,->] (2,0.) -- (5.8,0.);
\draw (2.2,0.5)node[right]{Solution of SDP~(\ref{eq:OptimalityProperty}),};
\draw (2.2,-0.5)node[right]{after standardisation,};
\draw (2.2,-1)node[right]{by using Algorithm~\ref{Alg1}.};
\draw (6.,1)node[right]{\bf Output 1};
\draw (6.,0.5)node[right]{We obtain $H_\Xi\in \mathbb{R}^{N\times r_0}$};
\draw (6.,0.)node[right]{s.t. $\rho^\star = H_\Xi H_\Xi^\intercal$.};
\draw (6.,-0.5)node[right]{svd: $H_\Xi =U\Sigma V^\intercal$ yields};
\draw (6.,-1)node[right]{the embedding $\Xi\in \mathbb{R}^N$.};
\draw (10.,0.5)node[right]{Extension};
\draw (10.,-0.5)node[right]{cfr.~(\ref{eq:ExtensionOnePoint}).};
\draw[thick,->] (10.,0.) -- (11.8,0);
\draw (11.8,1)node[right]{\bf Output 2};
\draw (11.8,0.5)node[right]{$\bar{\Xi}:X\to \mathbb{R}$ and};
\draw (11.8,0.)node[right]{$\bar{\rho}:X\times X\to \mathbb{R}$};
\draw (11.8,-0.5)node[right]{data-dependent};
\draw (11.8,-1)node[right]{p.s.d. kernel.};

\end{tikzpicture}
\caption{Schematic view of the numerical strategy for obtaining the SDP embedding and the data-dependent kernel~\cite{Wu2002}.\label{Fig:NumStrat}}
\end{figure}
It is worth mentioning that the method presented in this paper is more adapted to the case of homogeneous attributes. In the case of heterogeneous attributes, a substantial improvement is expected with a more appropriate choice of kernel or with the combination of a feature selection method. Incidentally, methods for dealing with heterogeneous features were proposed in the literature such as the Laplacian pyramids~\cite{Rabin}.

Empirically, the rank of the optimal solution $\rho^\star$ is often lower than $5$. Furthermore, the visualization obtained seems to be robust to the presence of data points of small degree (similar to outliers).
In order to illustrate the effect of points with small degrees, the data set similar to the clusters of Figure~\ref{Fig:Favorable} is now contaminated by outliers as shown in Figure~\ref{Fig:OutliersBandwidth} (top left), where the first largest eigenvalues of the diffusion kernel~(\ref{eq:DiffKernel}) are also displayed as a function of the bandwidth of the Gaussian kernel. We firstly observe that the spectral truncation is more reliable for a larger bandwidth since the spectral gap is more clear.  However, the outliers (stars) are not well separated from the clusters for such a large value of $\sigma$.
The effect of these outliers for different values of the bandwidth parameter is illustrated in Figure~\ref{Fig:OutliersBandwidth}, where the method proposed here and the diffusion map method are compared. In every case, the rank of the optimal solution $\rho^\star$ is equal to $2$ and the outliers remain for from the dense clusters as the bandwidth parameter varies. Namely, we observe qualitatively in the example of  Figure~\ref{Fig:OutliersBandwidth} that the result of the visualization does not significantly changes for different values of $\sigma$.
The embedding method is also illustrated on the Iris dataset in Figure~\ref{Fig:Iris} where the result is qualitatively similar to Diffusion Maps.\begin{figure}[H]
\setlength\figureheight{0.2\textwidth}
\setlength\figurewidth{0.3\textwidth}
\begin{minipage}{0.5\textwidth}
\input{SDPIrissigma1.tikz}
\end{minipage}
\begin{minipage}{0.5\textwidth}
\input{DMIrissigma1.tikz}
\end{minipage}
\caption{Comparison of the SDP embedding and the Diffusion Maps for the Iris dataset with $\sigma = 1$. The optimal solution of the SDP embedding satisfies $\rank(\rho^\star) = 2$. Each color refers a different class of points. \label{Fig:Iris}}
\end{figure} Here also the rank of the optimal solution is equal to $2$. In the case of the Wine dataset given in Figure~\ref{Fig:Wine}, we observe again the result of the SDP embedding gives still an interesting information (bottom left) when the bandwidth is chosen so small that Diffusion Maps only emphasizes the outliers (bottom right).
\begin{figure}[H]
\setlength\figureheight{0.2\textwidth}
\setlength\figurewidth{0.3\textwidth}
\begin{minipage}{0.5\textwidth}
\input{SDPwine10.tikz}
\end{minipage}
\begin{minipage}{0.5\textwidth}
\input{DMwine10.tikz}
\end{minipage}
\begin{minipage}{0.5\textwidth}
\input{SDPwine1c5.tikz}
\end{minipage}
\begin{minipage}{0.5\textwidth}
\input{DMwine1c5.tikz}
\end{minipage}
\caption{Comparison of the SDP embedding ($\rank(\rho^\star) = 2$) and the Diffusion Maps for the Wine dataset after a standardisation of the data points. Several different values of the bandwidth $\sigma$ are used.  Each color refers a different class of points.\label{Fig:Wine}}
\end{figure}

In Appendix, the reader will find additional illustrations of the SDP embedding method.

\subsection{Illustration of the out-of-sample extension}
The out-of-sample formula~(\ref{eq:Out-of-sample}) is illustrated in Figure~\ref{Fig:OOS} on the ``three clusters'' dataset of Figure~\ref{Fig:Favorable}. In order to have a reliable out-of-sample extension, the bandwidth of the Gaussian kernel should not be chosen too small.
Similarly, the digits $4$ and $5$ of the test set of the MNIST dataset are visualized in Figure~\ref{Fig:OOS_MNIST45} (left) thanks to the SDP embedding while the out-of-sample formula is used in order to embed the training set as displayed on the right of Figure~\ref{Fig:OOS_MNIST45}. Additional examples are given in~\ref{App:AddExamples}.
\begin{figure}[H]
\setlength\figureheight{0.2\textwidth}
\setlength\figurewidth{0.3\textwidth}
\begin{minipage}{0.5\textwidth}
\input{original3clusters.tikz}
\end{minipage}
\begin{minipage}{0.5\textwidth}
\input{embedded3clusters.tikz}
\end{minipage}
\caption{Three points are embedded thanks to the out-of-sample extension formula~(\ref{eq:Out-of-sample}) (black triangle, square and diamond). The color refers to the probability of presence at stationarity $\phi_0$. We found  $\rank(\rho^\star) = 2$. (The initial data is on the lhs and the embedding on the rhs.) \label{Fig:OOS}}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}{0.45\textwidth}
\centering\textbf{SDP Embedding}\\
\includegraphics[scale=0.6,trim={4cm 8.5cm 4cm 9.cm},clip]{MNIST45_Test_bw4.pdf}
\end{minipage}
\begin{minipage}{0.45\textwidth}
\centering\textbf{Out-of-sample extension}\\
\includegraphics[scale=0.6,trim={4cm 8.5cm 4cm 9.cm},clip]{MNIST45_OOSTrain_bw4.pdf}
\end{minipage}
\caption{On the left, SDP embedding for the ``test'' part of the MNIST dataset  for the digits $4$ and $5$ including in total 1874 points (left, blue and red stars) with a bandwidth $\sigma =4$. We found $\rank(\rho^\star) = 2$. On the right, the ``training'' part of the MNIST dataset (11263 points, blue and red dots) is embedded thanks to the extension formula~(\ref{eq:Out-of-sample}), while the stars are the 1874 points of the test set. Concerning the terminology, notice that the smallest dataset (``test'' part) is used for the training problem, i.e., solving the SDP.
\label{Fig:OOS_MNIST45}}
\end{figure}

\FloatBarrier
\section{Conclusion}

We may conclude as follows. We can start with a kernel $K:X\times X\to \mathbb{R}$ satisfying Assumption~\ref{Assumption1}, given for instance by
\[
K(x,y) = \frac{k(x,y)}{\sqrt{\bar{d}_\Omega(x)}\sqrt{\bar{d}_\Omega(y)}}-\sqrt{\frac{\bar{d}_\Omega(x)}{\vol(\Omega)}}\sqrt{\frac{\bar{d}_\Omega(y)}{\vol(\Omega)}},
\]
which is in general not positive semi-definite on $X\times X$ but only when restricted to $\Omega\times \Omega$, then
we obtain by solving~(\ref{eq:OptimalityProperty}) a Mercer kernel $\rho:\Omega\times \Omega\to \mathbb{R}$ of low rank which can be extended to a Mercer kernel on $X\times X$. The numerical method yields in fact the eigendecomposition $\rho = \Xi \Xi^\intercal$ which can be used in order to embed the data points in $\Omega$. A natural out-of-sample formula is also provided. It is observed numerically that the embedding is robust to the presence of outliers. The schematic numerical strategy is summarized in Figure~\ref{Fig:NumStrat}.\\

{\bf Acknowledgements}\\
\noindent {\small
M.F. acknowledges stimulating discussions with A. Themelis.
The authors also thank the following organizations. EU: The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC AdG A-DATADRIVE-B (290923). This paper reflects only the authors' views, the Union is not liable for any use that may be made of the contained information.
Research Council KUL: GOA/10/09 MaNet, CoE PFV/10/002 (OPTEC),
BIL12/11T; PhD/Postdoc grants
Flemish Government:
FWO: projects: G.0377.12 (Structured systems), G.088114N (Tensor based
data similarity); PhD/Postdoc grants
IWT: projects: SBO POM (100031); PhD/Postdoc grants
iMinds Medical Information Technologies SBO 2014
Belgian Federal Science Policy Office: IUAP P7/19 (DYSCO, Dynamical
systems, control and optimization, 2012-2017). This paper presents research results of the Concerted Research Action (ARC) programme supported by the Federation Wallonia-Brussels (contract ARC 14/19-060 on Mining and Optimization of Big Data Models)}\\

redpajama_set_name: RedPajamaArXiv