#  Online landmark replacement for out-of-sample dimensionality reduction
methods

Chanon Thongprayoon1, Naoki Masuda1,2,3,∗

###### Abstract

A strategy to assist visualization and analysis of large and complex data sets
is dimensionality reduction, with which one maps each data point into a low-
dimensional manifold. However, various dimensionality reduction techniques are
computationally infeasible for large data. Out-of-sample techniques aim to
resolve this difficulty; they only apply the dimensionality reduction
technique on a small portion of data, referred to as landmarks, and determine
the embedding coordinates of the other points using landmarks as references.
Out-of-sample techniques have been applied to online settings, or when data
arrive as time series. However, existing online out-of-sample techniques use
either all the previous data points as landmarks or the fixed set of landmarks
and therefore are potentially not good at capturing the geometry of the entire
data set when the time series is non-stationary. To address this problem, we
propose an online landmark replacement algorithm for out-of-sample techniques
using geometric graphs and the minimal dominating set on them. We
mathematically analyze some properties of the proposed algorithm, particularly
focusing on the case of landmark multidimensional scaling as the out-of-sample
technique, and test its performance on synthetic and empirical time series
data.

1Department of Mathematics, State University of New York at Buffalo, Buffalo,
NY, USA.
2Institute for Artificial Intelligence and Data Science, State University of
New York at Buffalo, Buffalo, NY, USA
3Center for Computational Social Science, Kobe University, Kobe, Japan
*Corresponding author<EMAIL_ADDRESS>
##  1 Introduction

One way to aid visualization and analysis of high-dimensional data is to seek
representations of the given set of data points in low-dimensional spaces,
including the case of non-Euclidean manifolds [23, 15, 66]. Finding a low-
dimensional representation of high-dimensional data is referred to as
dimensionality reduction [33, 28, 47]. Examples of dimensionality reduction
methods are the principal component analysis (PCA) [34], multidimensional
scaling (MDS) [22], locally linear embedding (LLE) [10], isometric feature
mapping (ISOMAP) [11], and $t$-distributed stochastic neighbor embedding
($t$-SNE) [29]. These and many other methods aim to preserve the dissimilarity
among data points as much as possible through the dimensionality reduction
process [33, 48]. Applications of dimensionality reduction include pattern
recognition, statistical analysis of multi-variable data, complexity
reduction, feature extraction, and signal processing, to name a few [23, 4,
65, 76].

In general, running a dimensionality reduction algorithm may require a
prohibitive amount of time when a data set is large. For example, the PCA,
MDS, and ISOMAP have the cubic time complexity [33], while other methods such
as LLE and $t$-SNE have the quadratic time complexity [33, 29]. One strategy
for suppressing the time complexity of dimensionality reduction algorithms is
to approximate the embedding coordinates of the data points. To do so, one can
turn to approximate minimization of the cost function of an embedding method
[52, 40, 43]. Alternatively, one can construct a computationally feasible
projection matrix, whose entries are random variables that take values from
probability distributions of data points of one’s choice, such as the Gaussian
distribution [9, 37, 35] or the distributions introduced in [17, 12]. Such a
projection matrix applied to the data points in the original high-dimensional
Euclidean space provides their approximate embedding coordinates in a low-
dimensional Euclidean space [17, 52, 35, 12, 9, 37]. Another strategy is to
exactly embed only a relatively small subset of the entire data, referred to
as out-of-sample techniques [50, 19, 74, 63, 49, 33]. An out-of-sample
dimensionality reduction method calculates the precise embedding coordinate
only for a small portion of data points, called the landmarks, and
approximates the coordinates for the remaining data points. Because one only
approximately calculates the coordinates of the non-landmark data points, the
embedding of the entire data set with an out-of-sample method usually runs
faster than the original embedding method. Examples of out-of-sample
techniques are the landmark multidimensional scaling (LMDS) [15, 20], landmark
isometric feature mapping (L-ISOMAP) [15], landmark diffusion maps [63],
kernel $t$-SNE [50], and bi-kernel $t$-SNE [74].

We have discussed embeddings of static data sets. In practice, one may want to
embed a stream of data arriving in real time into a low-dimensional space.
Out-of-sample techniques are directly applicable to such online settings
because one can calculate the embedding coordinates of newly arriving data
points, which are regarded as out-of-sample data. In general, online methods
should suppress the computational time by avoiding repeatedly applying a
dimensionality reduction procedure every time a new data point arrives [72,
71]. Known online dimensionality reduction techniques include online PCA [72,
7, 14, 59, 46], online Laplacian eigenmaps [55], online ISOMAP [21, 25],
online L-ISOMAP [25], and online LLE [30]. There are also downstream tasks
facilitated by online dimensionality reduction techniques for stream data,
such as online anomaly detection for cybersecurity [51], online classification
and clustering of data points [56, 42, 54, 44], and autonomous control in
robotics [53].

The online dimensionality reduction techniques mentioned above are instances
of out-of-sample techniques. Regarding the selection of landmarks, there are
two major classes of these methods. First, some algorithms use all the
previous data points as landmarks and therefore the number of landmarks keeps
growing [7, 14, 59, 46, 72, 21, 25, 55, 30]. This class of methods does not
restrict the number of landmarks and hence does not particularly aim to reduce
runtime. Second, other algorithms use a set of landmarks that is fixed over
time once it is initialized [20, 50, 63, 74, 25]. An earlier study suggested
that using evenly distributed landmarks preserves the data structure better
than using a set of poorly distributed landmarks [20]. Therefore, one idea is
to choose landmarks from the data set uniformly at random. However, the new
data points may explore an unexplored part of the space of the data, and they
may be far away from any of the already revealed data points. In this case,
embedding with an initially fixed set of landmarks may cause substantial
distortion because they are unevenly distributed given the data set including
the newer data points. It may be then better to adaptively replace some
landmarks as new data points arrive to assure that all landmarks do not
congregate in a small region of the data space and the set of landmarks
provides a reasonable scaffold of the entire data set at any given time. We
currently lack such landmark replacement algorithms for online out-of-sample
embedding.

We remark that there are online dimensionality reduction techniques that do
not use out-of-sample techniques. In fact, they require one to increase the
dimension of the embedding space as new data points arrive [72, 7, 14, 46],
which may be inconvenient in practice.

To fill this gap, we study online landmark selection for out-of-sample
techniques. By constructing an adaptive geometric graph [18, 36] from the
given data set and exploiting an algorithm that finds the online dominating
set [62], we propose an algorithm of online landmark replacement for out-of-
sample dimensionality reduction methods. Upon the arrival of a new data point,
the proposed algorithm determines whether or not a current landmark should be
replaced. We then provide mathematical underpinnings of the proposed
algorithm, specifically in the case of the LMDS, and apply the proposed
algorithm to three multidimensional time series data sets.

##  2 Online landmark replacement algorithm

In this section, we propose an online landmark replacement algorithm for out-
of-sample embedding methods. We introduce building blocks for our algorithm in
sections 2.1 and 2.2. Then, we explain our online algorithm in section 2.3.

###  2.1 Geometric graph, dominating set, and landmark

Figure 1: A graph with 12 nodes. The set $L=\\{v_{3},v_{6},v_{10}\\}$, shown
in blue, is a dominating set because, for each $r\in\\{1,\ldots,12\\}$,
$v_{r}$ either belongs to $L$ or is adjacent to an element in $L$. On the
other hand, $L_{0}=\\{v_{4},v_{7},v_{11}\\}$, shown in red, is not a
dominating set because $v_{1},v_{9},$ and $v_{12}$ do not belong to $L_{0}$ or
is adjacent to an element in $L_{0}$. Here,
$OnlyBy(v_{3})=\\{v_{1},v_{2},v_{3},v_{4}\\}$. For instance, $v_{5}\not\in
OnlyBy(v_{3})$ because $v_{5}$ is dominated by both $v_{3}$ and $v_{6}$.

In this section, we introduce the notion of landmark via the definition of
geometric graph and dominating set. We let $V=\\{x_{1},x_{2},\ldots,x_{n}\\}$
be the set of data, arriving at $0\leq t_{1}<t_{2}<\cdots<t_{n}$, $d:V\times
V\rightarrow\mathbb{R}$ be a metric, and $\varepsilon$ be a positive number.
We define a network with node set $V$ and edge set given by
$E=\\{(x_{r},x_{s}):d(x_{r},x_{s})<\varepsilon\\}$. This network is called a
geometric graph [18, 36]. We also denote by $L$ $(\subset V)$ the set of
landmarks. Our problem of landmark replacement is to update $L$ that contains
at most $m$ landmarks given $V$ and $d(x_{r},x_{s})$ $\forall
r,s\in\\{1,\ldots,n\\}$ such that $L$ has some desirable properties. We fix
$m$.

To find $L$, we turn to the dominating set problem [58, 2]. The set
$L\subseteq V$ is said to be a dominating set of a graph $(V,E)$ if, for each
node $v$, either $v$ belongs to $L$ or $v$ is adjacent to a node in $L$. In
the former case, we say that $v$ covers itself (see Fig. 1 for an
illustration). In the latter case, we say that $v$ is covered by a different
node in $L$. If $v$ neither belongs to $L$ nor is covered by a different node
in $L$, then we say that $v$ is uncovered. If there is at least one $v$ that
is uncovered, then $L$ is not a dominating set. We aim to select a dominating
set $L$ to be used as the set of landmarks for out-of-sample embedding
methods.

###  2.2 Selecting landmarks offline

Here, we introduce a procedure to determine $L$ in an offline setting. We
assume that $V$ is given and fixed. If $n\leq m$, then we set $L=V$.
Otherwise, we apply Algorithm 1 to find a dominating set $L\subseteq V$ for
network $(V,E)$ [3] and use $L$ as the set of landmarks.

Algorithm 1 Greedy dominating set

1:Set $L=\emptyset$

2:while there is an uncovered node do

3: Choose the $x_{r}\notin L$ such that $x_{r}$ covers the largest number of
uncovered nodes.

4: $L\leftarrow L\cup\\{x_{r}\\}$

5:end while

6:return $L$

###  2.3 Updating landmarks online

We now aim to update the set of landmarks, $L$, assuming that the data points
sequentially arrive. When a new data point $x_{r}$ arrives, we either keep the
current $L$ or replace an element in $L$ by another element.

There are algorithms that allow updating of the dominating set as nodes are
sequentially added or deleted [16, 70, 5, 62]. For instance, one method
recomputes a dominating set on time-varying unit disk graphs upon each
addition of deletion of a node [16]. In this algorithm, there is no limitation
to the size of the dominating set, which we denote by $\left|L\right|$. If the
added node is covered by the current dominating set, then no action is
required. Otherwise, we include the added node into the dominating set. When
deleting a node $v$ from the network, if $v$ does not belong to the dominating
set, no action is taken. Otherwise, there may be uncovered nodes as a result
of deleting $v$ from the dominating set. If the latter is the case, we add all
the uncovered neighbors of $v$ to the dominating set. This last step may
drastically increase $\left|L\right|$. Another online algorithm, which is the
latter of the two algorithms proposed in [62], seeks to find a minimal
dominating set and does not restrict $m$. In addition to allowing sequentially
adding nodes to $L$, this algorithm permits sequential addition of edges to
the network, which is a useful property for our landmark replacemenet
algorithm as we explain later. Two other online algorithms do not allow us to
add an edge between the existing nodes [5]. Reference [70] provides lower and
upper bounds on the size of dominating sets for deterministic and randomized
online dominating set algorithms.

All these algorithms allow $\left|L\right|$ to indefinitely grow and do not
provide a mechanism to control $\left|L\right|$. Because we want to maintain
the dominating set to satisfy $\left|L\right|\leq m$ as we add nodes, we
propose to allow $\varepsilon$ to adaptively increase, which is equivalent to
adding edges to the network, with the aim of decreasing $\left|L\right|$ when
$\left|L\right|$ exceeds $m$. Although the two algorithms proposed in [5],
discussed above, allow us to add edges, one can only add new edges incident to
any new node $v$ when $v$ has joined the network. These algorithms do not
allow us to add an edge $(v,u)$ at a later stage, which we may want to do in
our landmark replacement procedure in which we increase $\varepsilon$ as we
add nodes, as we describe below. Therefore, we need to resort to a different
online dominating set algorithm to realize our goal.

We adapt an online algorithm that finds the minimal dominating set of time-
varying networks [62] to our landmark replacement problem. With our proposed
method, we control the number of landmarks in a dynamic manner by gradually
increasing $\varepsilon$, which corresponds to adding edges to the geometric
graph.

We assume that $V=\\{x_{1},\ldots,x_{n}\\}$ and that $x_{n+1}$ is a new data
point. We regard that $x_{n+1}$ is a new node of the network, which directly
connects to each of the existing nodes within distance $\varepsilon$ from
$x_{n+1}$ by a new edge. Then, exactly one of the following three statements
holds true:

  1. 1.

$\mathcal{N}(x_{n+1})\cup\\{x_{n+1}\\}\neq\\{x_{n+1}\\}$ and
$\left(\mathcal{N}(x_{n+1})\cup\\{x_{n+1}\\}\right)\cap L\neq\emptyset$,

  2. 2.

$\mathcal{N}(x_{n+1})\cup\\{x_{n+1}\\}\neq\\{x_{n+1}\\}$ and
$\left(\mathcal{N}(x_{n+1})\cup\\{x_{n+1}\\}\right)\cap L=\emptyset$,

  3. 3.

$\mathcal{N}(x_{n+1})\cup\\{x_{n+1}\\}=\\{x_{n+1}\\}$,

where $\mathcal{N}(x_{n+1})$ is the set of neighbors of $x_{n+1}$.

In the first case, we do not need to update $L$ given the current value of
$\varepsilon$. Therefore, there is no need to update $\varepsilon$. To discuss
the other two cases, we define $OnlyBy(x)$ for any $x\in L$ as follows:

$OnlyBy(x)=\\{y\in\mathcal{N}(x)\cup\\{x\\}:\left|\left(\mathcal{N}(y)\cup\\{y\\}\right)\cap
L\right|=1\\}.$ (1)

In other words, $OnlyBy(x)$ is the set of neighbors of $x$ that only $x$
dominates (see Fig. 1 for an illustration).

Let us consider the second case. If $\left|L\right|<m$, then we add $x_{n+1}$
to the set of landmarks, $L$. Otherwise, because $\left|L\right|=m$, we use
Algorithm 2 to determine if we replace a landmark in $L$ by $x_{n+1}$ or by
another non-landmark data point. Reference [62] provides an algorithm for
adding and removing nodes and/or edges, including a criterion to decide
whether or not we add $x_{n+1}$ to $L$. Unlike the original work, we limit the
size of $L$ by imposing $\left|L\right|\leq m$, and we cannot add an arbitrary
edge to the network because the edges to be added depends on the proximity
between a pair of nodes measured by $d$. Instead, we gradually increase the
$\varepsilon$ value to add an edge one by one in Algorithm 2. In summary,
Algorithm 2 either adds $x_{n+1}$ or a neighbor of $x_{n+1}$ to $L$; we make
the decision based on $\deg(x_{n+1})$, which is the degree of node $x_{n+1}$.
Then, we add an edge between a pair of nodes that realizes the smallest
distance among all the node pairs with distance greater than the current
$\varepsilon$. This entails an increase in $\varepsilon$. If neither of the
two nodes joined by the new edge is an element of $L$, then we further
increase $\varepsilon$ such that the next pair of nodes is additionally joined
by an edge. We repeat adding edges in this manner until either of the two
nodes that are joined by the new edge is a landmark. Finally, to decide which
landmark will be replaced, we call Algorithm 3 (i.e., function REMOVE) [62].
Algorithm 3 returns a landmark after the removal of which $L$ remains a
dominating set. If such a landmark exists, then Algorithm 2 removes this
landmark. Otherwise, i.e., if Algorithm 3 returns an empty set, Algorithm 2
sequentially keeps adding edges by gradually increasing $\varepsilon$ until
REMOVE returns a landmark that can be removed. Lines 7-12 and 14-18 in
algorithm 2 are novel, and the remainder of Algorithm 2 including Algorithm 3
is the same as the original one [62].

In the third case, $x_{n+1}$ is an isolated node in $(V,E)$. This is similar
to the second case in that $x_{n+1}$ is not adjacent to any element in $L$. If
$\left|L\right|<m$, then we just add $x_{n+1}$ to $L$. Otherwise, i.e., if
$\left|L\right|=m$, we apply Algorithm 2 to $(V,E)$. Because
$\deg(x_{n+1})=0$, lines 1 and 2 of Algorithm 2 imply that one adds $x_{n+1}$
to $L$. Then, because there are now $m+1$ landmarks, we remove one of them by
carrying out lines 7-18 in Algorithm 2.

Algorithm 2 Landmark replacement

Input $(V,E)$

1:if $\deg(x_{n+1})\leq 2\sqrt{\left|E\right|}$ then

2: $L\leftarrow L\cup\\{x_{n+1}\\}$

3:else

4: Choose an arbitrary $\tilde{x}\in\mathcal{N}(x_{n+1})$ such that
$\deg(\tilde{x})\leq\sqrt{\left|E\right|}$

5: $L\leftarrow L\cup\\{\tilde{x}\\}$

6:end if

7:$\left(x_{r},x_{s}\right)=\text{argmin}_{\\{x_{r^{\prime}},x_{s^{\prime}}\\}\not\in
E}d(x_{r^{\prime}},x_{s^{\prime}})$

8:Update $\varepsilon=d(x_{r},x_{s})$ and $E\leftarrow
E\cup\\{(x_{r},x_{s})\\}$

9:while $x_{r}\not\in L\text{ and }x_{s}\not\in L$ do

10:
$\left(x_{r},x_{s}\right)=\text{argmin}_{\\{x_{r^{\prime}},x_{s^{\prime}}\\}\not\in
E}d(x_{r^{\prime}},x_{s^{\prime}})$

11: Update $\varepsilon=d(x_{r},x_{s})$ and $E\leftarrow
E\cup\\{(x_{r},x_{s})\\}$

12:end while

13:$\pi=\text{REMOVE}(x_{r},x_{s})$

14:while $\pi=\emptyset$ do

15: Update $\varepsilon=\min_{\\{x_{r^{\prime}},x_{s^{\prime}}\\}\not\in
E}d(x_{r^{\prime}},x_{s^{\prime}})$ and $E\leftarrow
E\cup\\{(x_{r^{\prime}},x_{s^{\prime}})\\}$

16: $\pi=\text{REMOVE}(x_{r},x_{s})$

17:end while

18:$L\leftarrow L\setminus\pi$

19:return $L,\varepsilon$

Algorithm 3 REMOVE

Input $x_{r},x_{s}$

1:if $x_{r}\in L$ then

2: if $OnlyBy(x_{r})=\emptyset$ then

3: return $x_{r}$

4: else if $x_{s}\in L\text{ and }OnlyBy(x_{s})=\emptyset$ then

5: return $\\{x_{s}\\}$

6: else

7: return $\emptyset$

8: end if

9:else

10: if $OnlyBy(x_{s})=\emptyset$ then

11: return $\\{x_{s}\\}$

12: else

13: return $\emptyset$

14: end if

15:end if

Algorithm 4 Two-way merge sorting algorithm [8]

Input $X=\left(X[1],\ldots,X[p]\right),Y=\left(Y[1],\ldots,Y[q]\right)$, where
$X[r^{\prime}]\leq X[r^{\prime}+1]$ $\forall r^{\prime}\in\\{1,\ldots,p-1\\}$,
$Y[s^{\prime}]\leq Y[s^{\prime}+1]$ $\forall s^{\prime}\in\\{1,\ldots,q-1\\}$

1:Initialize an empty array $Z=()$.

2:$r=1,s=1$

3:while $r<p+1\text{ and }s<q+1$ do

4: if $X[r]\leq Y[s]$ then

5: $Z.\text{append}\left(X[r]\right)$

6: $r=r+1$

7: else

8: $Z.\text{append}\left(Y[s]\right)$

9: $s=s+1$

10: end if

11:end while

12:if $r=p+1$ then

13: return $\text{Concatenate}\left(Z,\left(Y[s],\ldots,Y[q]\right)\right)$

14:else

15: return $\text{Concatenate}\left(Z,\left(X[r],\ldots,X[p]\right)\right)$

16:end if

##  3 Mathematical properties

In this section, we show some mathematical properties related to the proposed
landmark replace algorithm, in particular its time complexity.

######  Theorem 3.1.

[3] For an undirected graph $(V,E)$,

$\left|L\right|\leq\left|V\right|+1-\sqrt{2\left|E\right|+1},$ (2)

where $L$ is the dominating set of $(V,E)$ obtained from Algorithm 1.

The following Corollary of Theorem 3.1 tells us how we should set the value of
$\varepsilon$ to guarantee that the number of landmarks obtained by Algorithm
1 is at most $m$.

######  Corollary 3.2.

Consider an undirected graph $(V,E)$, where $V=\\{x_{1},\ldots,x_{n}\\}$. If
we set the value of parameter $\varepsilon$ such that

$\left|E\right|=\left\lceil\frac{(n-m+1)^{2}-1}{2}\right\rceil,$ (3)

where $E=\\{(x_{r},x_{s}):d(x_{r},x_{s})<\varepsilon\\}$, then $L$ obtained by
Algorithm 1 satisfies $\left|L\right|\leq m$.

###### Proof.

Using Eq. (3), we obtain

$\displaystyle\sqrt{2\left|E\right|+1}$
$\displaystyle=\sqrt{2\left(\left\lceil\frac{(n-m+1)^{2}-1}{2}\right\rceil\right)+1}$
$\displaystyle\geq\sqrt{2\left(\frac{(n-m+1)^{2}-1}{2}\right)+1}$
$\displaystyle=n-m+1.$ (4)

Equation (4) implies that

$n+1-\sqrt{2\left|E\right|+1}\leq m.$ (5)

By combining Theorem 3.1 and Eq. (5), we obtain

$\left|L\right|\leq n+1-\sqrt{2\left|E\right|+1}\leq m.$ (6)

∎

######  Theorem 3.3.

Let $x$ be a node in an undirected graph $(V,E)$. If
$\deg(x)>2\sqrt{\left|E\right|}$, then there is a neighbor $x^{\prime}$ of $x$
such that $\deg(x^{\prime})\leq\sqrt{\left|E\right|}$.

###### Proof.

Suppose to the contrary that $\deg(x^{\prime})>\sqrt{\left|E\right|}$ for all
$x^{\prime}\in\mathcal{N}(x)$. Therefore, we obtain

$\displaystyle 2\left|E\right|$
$\displaystyle\geq\sum_{x^{\prime}\in\mathcal{N}(x)}\deg(x^{\prime})$
$\displaystyle>\sum_{x^{\prime}\in\mathcal{N}(x)}\sqrt{\left|E\right|}$
$\displaystyle=\sqrt{\left|E\right|}\cdot\left|\mathcal{N}(x)\right|$
$\displaystyle>2\left|E\right|,$ (7)

which is a contradiction. ∎

######  Remark 3.4.

This theorem extends a theorem given in [62], which provides a sufficient
condition $\deg(x)>2\sqrt{\left|E\right|}+1$ rather than
$\deg(x)>2\sqrt{\left|E\right|}$.

######  Remark 3.5.

A proof of the theorem in [62], which is omitted in [62], is similar to this
proof.

Next, we formally define the Big $O$ notation to discuss time complexity of
the proposed algorithm in this section and other properties in section 5.

######  Definition 1.

Let $A$ be an unbounded subset of $\mathbb{R}^{+}$ such that
$f_{i}:A\rightarrow\mathbb{R}$ (with $i\in\\{1,2\\}$) is continuous. Then,
$f_{1}(x)=O\left(f_{2}(x)\right)$ as $x\rightarrow\infty$ if there exist
$M,C>0$ such that $\left|f_{1}(x)\right|\leq M\left|f_{2}(x)\right|$, $\forall
x>C$.

Alternatively, let $B$ be a subset of $\mathbb{R}$ such that
$g_{i}:B\rightarrow\mathbb{R}$ (with $i\in\\{1,2\\}$) is continuous.
Furthermore, let $x_{0}$ be a real number. Then, $g_{1}(x)=O(g_{2}(x))$ as
$x\rightarrow x_{0}$ if there exist $\delta,M>0$ such that
$\left|g_{1}(x)\right|\leq M\left|g_{2}(x)\right|$, $\forall x$ satisfying
$0<\left|x-x_{0}\right|<\delta$.

######  Lemma 3.6.

Assume that we have constructed the adjacency list of graph $(V,E)$, where
$V=\\{x_{1},\ldots,x_{n}\\}$, and
$E=\\{(x_{r},x_{s}):d(x_{r},x_{s})<\epsilon\\}$, where
$r,s\in\\{1,\ldots,n\\}$; the adjacency list is the list of neighbors of each
node. Then, the time complexity for computing $OnlyBy(z)$ for a given $z\in L$
is $O(mn)$.

###### Proof.

First, if $z$ is an isolated node, then $\mathcal{N}(z)\cup\\{z\\}=\\{z\\}$
and $\left|\left(\mathcal{N}(z)\cup\\{z\\}\right)\cap
L\right|=\left|\\{z\\}\right|=1$. Therefore, we obtain $OnlyBy(z)=\\{z\\}$,
and the time complexity for computing $OnlyBy(z)$ is $O(1)$. Now, we assume
that $z$ is not isolated. Then, we check if $\mathcal{N}(z)\subseteq L$. To do
so, we examine the list of neighbors of $z$. If $\mathcal{N}(z)\subseteq L$,
then every neighbor of $z$ is a landmark, and we obtain $OnlyBy(z)=\emptyset$.
The time complexity for checking whether or not $\mathcal{N}(z)\subseteq L$ is

$O\left(\left|L\right|\deg(z)\right)\leq O\left(mn\right).$ (8)

If $\mathcal{N}(z)\nsubseteq L$, then we create two arrays of length $n$,
denoted by $A=(A[1],\ldots,A[n])$ and $B=(B[1],\ldots,B[n])$ and initialize
each entry of $A$ and $B$ by zero. Then, we set $A[r]=1$ for each
$x_{r}\in\mathcal{N}(z)\cup\\{z\\}$. Similarly, we set $B[r]=1$ for each
$x_{r}$ satisfying $x_{r}\in\mathcal{N}(z^{\prime})\cup\\{z^{\prime}\\}$ for
any $z^{\prime}\in L\setminus z$. The time complexity of setting $A$ and $B$
is $O(n)$ and $O(mn)$, respectively. Because every node is covered by a
landmark, either $A\left[r\right]=1$ or $B\left[r\right]=1$ holds true for all
$r\in\\{1,\ldots,n\\}$. Therefore, there are three cases to consider:
$(A\left[r\right],B\left[r\right])=(1,1)$, $(1,0)$, and $(0,1)$.

If $(A\left[r\right],B\left[r\right])=(1,1)$, then node $x_{r}$ is covered by
both $z$ and at least one node in $L\setminus\\{z\\}$. In this case, we do not
add $x_{r}$ to $OnlyBy(z)$. If $(A\left[r\right],B\left[r\right])=(1,0)$, then
we add $x_{r}$ to $OnlyBy(z)$ because $z$ is the only landmark that covers
$x_{r}$. If $(A\left[r\right],B\left[r\right])=(0,1)$, then we do not add
$x_{r}$ to $OnlyBy(z)$ because $x_{r}$ is not a neighbor of $z$ while $x_{r}$
is covered by at least one node in $L\setminus\\{z\\}$. Therefore, the total
time complexity for constructing $OnlyBy(z)$ from the arrays $A$ and $B$ is
$O(n)$.

By combining these time complexities, we find that the computation of
$OnlyBy\left(z\right)$ requires
$O(1)+O\left(mn\right)+O\left(n\right)+O\left(mn\right)+O\left(n\right)=O\left(mn\right)$
time. ∎

We now analyze the time complexity for executing Algorithm 2 upon the addition
of a single new data point, $x_{n+1}$. We assume that the pairwise distance
between the existing points has been calculated and sorted and that the
pairwise distance between each existing point and $x_{n+1}$ has been
calculated but not sorted yet. We proceed by sorting the array of pairwise
distance between each existing point and $x_{n+1}$. Then, we apply Algorithm 4
to the two sorted arrays to obtain the sorted array of the pairwise distances
of all points including $x_{n+1}$. Algorithm 4, known as the two-way merge
sorting algorithm, runs in $O(p+q)$ time, where $p$ and $q$ are the lengths of
each input array [8].

######  Theorem 3.7.

Upon the arrival of $x_{n+1}$, we assume that $d(x_{1},x_{n+1})$, $\ldots$,
$d(x_{n},x_{n+1})$ have been computed, that
$\\{d(x_{1},x_{2}),d(x_{1},x_{3}),\ldots,d(x_{n-1},x_{n})\\}$ has been sorted,
and that we have updated the adjacency list. Then, the time complexity of
executing Algorithm 2 with $V=\\{x_{1},\ldots,x_{n+1}\\}$ and
$E=\\{(x_{r},x_{s}):r,s\in\\{1,\ldots,n+1\\},d(x_{r},x_{s})<\varepsilon\\}$ as
input is $O\left(n^{2}+e_{n+1}mn\right)$, where $e_{n+1}$ is the number of
edges added upon executing Algorithm 2, and $m=\left|L\right|$, i.e., the
number of landmarks.

###### Proof.

First, we sort the array $\\{d(x_{1},x_{n+1}),\ldots,d(x_{n},x_{n+1})\\}$,
which requires $O(n\log(n))$ time. Then, we use Algorithm 4 to merge sort the
sorted array of $\\{d(x_{1},x_{n+1}),\ldots,d(x_{n},x_{n+1})\\}$ and the
sorted array of $\\{d(x_{r},x_{s}):r,s\in\\{1,\ldots,n\\},r\neq s\\}$. This
merge sort requires $O(n^{2})$ time because the length of the second array is
$O(n^{2})$. To summarize these two steps, we have sorted $d(x_{r},x_{s})$ with
$\forall r,\forall s\in\\{1,\ldots,n+1\\}$, $r\neq s$ in $C_{1}=O(n\log
n)+O(n^{2})=O(n^{2})$ time.

Theorem 3.3 implies that one only needs to inspect at most
$2\left\lceil\sqrt{\left|E\right|}\right\rceil+1$ neighbors of $x_{n+1}$ to be
able to obtain a neighbor whose degree is less than or equal to
$\sqrt{\left|E\right|}$. Therefore, we can carry out line 4 of Algorithm 2 in
$O\left(\sqrt{\left|E\right|}\right)$ time. Because we have computed and
sorted the pairwise distances between the nodes, lines 7–8 run in
$O\left(1\right)$ time. Therefore, the time complexity of carrying out lines
4–8 is $C_{2}=O\left(\sqrt{\left|E\right|}\right)$. We denote by
$e^{(1)}_{n+1}$ the number of edges that one adds to the network in lines
9–12. Moreover, we update the adjacency list upon adding $e^{(1)}_{n+1}$ edges
to the network. Therefore, the total cost of carrying out lines 9–12 is
$C_{3}=O\left(e^{(1)}_{n+1}\right)+O\left(2e^{(1)}_{n+1}\right)=O\left(e^{(1)}_{n+1}\right)$,
where $O\left(2e^{(1)}_{n+1}\right)$ comes from the time complexity for
updating the adjacency list. Lemma 3.6 guarantees that line 13 runs in
$C_{4}=O\left(m(n+1)\right)$ time.

In lines 14–17, for each added edge, we update the adjacency list and then
call Algorithm 3. Algorithm 3 requires that we (i) check whether or not one or
both input nodes belong to $L$ and then (ii) compute an $OnlyBy$. The cost of
(i) is $O\left(m\right)$. Lemma 3.6 indicates that the cost of (ii) is
$O(m(n+1))$. The total time complexity for updating the adjacency list and
carrying out lines 14–17 is
$C_{5}=O\left(2e^{(2)}_{n+1}\right)+O\left(e^{(2)}_{n+1}m\right)+O\left(e^{(2)}_{n+1}m(n+1)\right)=O\left(e^{(2)}_{n+1}m(n+1)\right)$,
where $e^{(2)}_{n+1}$ is the number of edges added upon the execution of lines
14–17.

Finally, we add $C_{1}$, $C_{2}$, $C_{3}$, $C_{4}$, and $C_{5}$, and use
$\left|E\right|\leq n(n-1)/2$ to find the total time complexity of executing
Algorithm 2 as follows:

$\displaystyle
O\left(n^{2}+\sqrt{\left|E\right|}+e^{(1)}_{n+1}+m(n+1)+e^{(2)}_{n+1}m(n+1)\right)$
$\displaystyle=O\left(n^{2}+e^{(1)}_{n+1}+e^{(2)}_{n+1}m(n+1)\right)$
$\displaystyle\leq
O\left(n^{2}+\left[e^{(1)}_{n+1}+e^{(2)}_{n+1}\right]m(n+1)\right)$
$\displaystyle=O\left(n^{2}+e_{n+1}mn\right).$ (9)

∎

######  Corollary 3.8.

Suppose that we sequentially add $x_{n+1},x_{n+2},\ldots,x_{n^{\prime}}$ with
possible landmark replacements by iterating Algorithm 2 with
$\left|L\right|=m$. We further assume that the distance between the new data
point, $x_{\ell}$, and each existing point, i.e., $d(x_{r},x_{\ell})$,
$r\in\\{1,\ldots,\ell-1\\}$, has been computed upon the addition of $x_{\ell}$
and that $\\{d(x_{1},x_{2}),d(x_{1},x_{3}),\ldots,d(x_{\ell-2},d_{\ell-1})\\}$
has been sorted. Then, the time complexity for iterating Algorithm 2 over
$x_{n+1},x_{n+2},\ldots,x_{n^{\prime}}$ is $O\left(m{n^{\prime}}^{3}\right)$.

###### Proof.

For each new $x_{\ell}$ with $\ell\in\\{n+1,\ldots,n^{\prime}\\}$, Theorem 3.7
implies that running Algorithm 2 requires
$O\left((\ell-1)^{2}+e_{\ell}m(\ell-1)\right)=O\left(\ell^{2}+e_{\ell}m\ell\right)$
time. Therefore, the time complexity of iterating Algorithm 2 from $x_{n+1}$
to $x_{n^{\prime}}$ is

$\displaystyle\sum_{\ell=n+1}^{n^{\prime}}O\left(\ell^{2}+e_{\ell}m\ell\right)$
$\displaystyle\leq\sum_{\ell=n+1}^{n^{\prime}}O\left(\ell^{2}\right)+O\left(mn^{\prime}\right)\sum_{\ell=n+1}^{n^{\prime}}O\left(e_{\ell}\right)$
$\displaystyle\leq
O\left({n^{\prime}}^{3}\right)+O\left(mn^{\prime}\right)O\left(\binom{n^{\prime}}{2}\right)$
$\displaystyle=O\left(m{n^{\prime}}^{3}\right).$ (10)

∎

##  4 Landmark multidimensional scaling

The landmark replacement method proposed in section 2.3 works for any out-of-
sample dimensionality reduction method as long as the distance between each
pair of data points in the original space is defined. In the remainder of this
article, we use the LMDS to evaluate the proposed landmark replacement method
mathematically and numerically.

###  4.1 Algorithms

In this section, we briefly review the LMDS [20], which is an out-of-sample
extension of the classical MDS.

The classical MDS works as follows. We consider the set of $n$ data points to
be embedded, $S=\\{x_{1},\ldots,x_{n}\\}$, with the distance between two data
points being given by $d(x_{r},x_{s})$. The MDS is a map from $S$ to
$\mathbb{R}^{k}$. To construct the MDS, we denote the $n\times n$ squared
distance matrix by $\Delta=[\Delta_{rs}]$, where
$\Delta_{rs}=d^{2}(x_{r},x_{s})$. We then define the double mean-centered dot
product matrix by $D=-\frac{1}{2}H_{n}\Delta H_{n}$, where
$H_{n}=I_{n}-\frac{1}{n}J_{n}$, matrix $I_{n}$ is the $n\times n$ identity
matrix, and $J_{n}$ is the $n\times n$ matrix in which all the entries are
equal to 1. The coordinates of $x_{r}$ in $\mathbb{R}^{k}$ after the mapping
is the $r$th column of

$L_{k}=\begin{bmatrix}\sqrt{\lambda_{1}}v_{1}^{\top}\\\ \vdots\\\
\sqrt{\lambda_{k}}v_{k}^{\top}\end{bmatrix},$ (11)

which is a $k\times n$ matrix. In Eq. (11), $\lambda_{r}$ is the $r$th largest
eigenvalue of $D$, which we assume to be positive, $v_{r}$ is the right
eigenvector associated with eigenvalue $\lambda_{r}$, normalized in terms of
the $\ell^{2}$-norm, and ⊤ represents the transposition. The time complexity
of the MDS is $O(n^{3})$ owing to the eigenvalue and eigenvector computation
[49]. Note that the MDS preserves the pairwise distance between all pairs of
data points in the Euclidean embedding space if and only if matrix
$D=-\frac{1}{2}H_{n}\Delta H_{n}$ is positive semidefinite [6, 13].

The LMDS only applies the MDS to a subset of the entire data set called the
landmarks [20, 15, 24]. The LMDS then determines the mapping of other
arbitrary data points with the help of the landmarks. We denote by
$L=\\{x_{i_{1}},\ldots,x_{i_{m}}\\}\subseteq S$ the set of landmarks and apply
the MDS to any $x\in L$ using Eq. (11). The coordinate in $\mathbb{R}^{k}$ of
an arbitrary out-of-sample data point $x\notin L$ is given by

$\psi(x)=-\frac{1}{2}L^{\prime}_{k}(\delta_{x}-\delta_{\mu}),$ (12)

where

$L^{\prime}_{k}=\begin{bmatrix}\frac{1}{\sqrt{\lambda_{1}}}v_{1}^{\top}\\\
\vdots\\\
\frac{1}{\sqrt{\lambda_{k}}}v_{k}^{\top}\end{bmatrix}\in\mathbb{R}^{k\times
m},$ (13) $\delta_{x}=\begin{bmatrix}d^{2}(x,x_{i_{1}})\\\ \ \vdots\\\
d^{2}(x,x_{i_{m}})\end{bmatrix},$ (14)
$\delta_{\mu}=\frac{1}{m}\sum_{r=1}^{m}\begin{bmatrix}d^{2}(x_{i_{1}},x_{i_{r}})\\\
\ \vdots\\\ d^{2}(x_{i_{m}},x_{i_{r}})\end{bmatrix}.$ (15)

###  4.2 Goodness of embedding

To assess the quality of the LMDS, we use the following two goodness of fit
criteria. First, the normalized stress function [22, 23, 1] is defined by

$\sigma=\sqrt{\frac{\sum_{r=1}^{n}\sum_{s=1}^{r-1}\left[d(x_{s},x_{r})-\left|\psi(x_{s})-\psi(x_{r})\right|_{2}\right]^{2}}{\sum_{r=1}^{n}\sum_{s=1}^{r-1}d^{2}(x_{s},x_{r})}}\,.$
(16)

Note that $\left|\psi(x_{s})-\psi(x_{r})\right|_{2}$ is the Euclidean distance
between the two data points in the embedding space. One obtains a smaller
$\sigma$ if the MDS better preserves the distance between pairs of data points
in the original space. A guideline of acceptable $\sigma$ values is
$\sigma<0.15$ [22].

Second, we denote by $\sigma_{L}$ the normalized stress computed only on the
basis of the set of landmarks. A small $\sigma_{L}$ value indicates that the
MDS well preserves the distance structure for the set of $m$ landmarks.

##  5 Dislocation of the embedding coordinates upon landmark replacement

In this section, we mathematically evaluate how much a data point moves in the
embedding space when one replaces a landmark with another data point. We
provide relevant basic definitions and theorems on linear operators in
Appendix A.

######  Definition 2.

Given a pair of $a\times b$ real-valued matrices $P=[p_{rs}]$ and
$Q=[q_{rs}]$, the Frobenius inner product is defined by

$\left\langle P,Q\right\rangle_{F}=\sum_{r=1}^{a}\sum_{s=1}^{b}p_{rs}q_{rs}.$
(17)

Additionally, we define the induced Frobenius norm by

$\left\|P\right\|_{F}=\sqrt{\left\langle P,P\right\rangle_{F}}$ (18)

and the Frobenius distance by

$d_{F}(P,Q)=\left\|P-Q\right\|_{F}.$ (19)

######  Theorem 5.1.

(De Silva and Tenenbaum [20]) Let $\Delta$ be the squared-distance matrix of
$m$ landmarks, and let $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{m}$
be the eigenvalues of $D=-\frac{1}{2}H_{m}\Delta H_{m}$ with $\lambda_{k}>0$
and $\lambda_{k}>\lambda_{k+1}$, where $k$ is the dimension of the embedding
space, $\mathbb{R}^{k}$. Consider a perturbation of $\Delta$ given by
$\tilde{\Delta}=\Delta+\varepsilon\Phi+O(\varepsilon^{2})$, where
$|\varepsilon|\ll 1$. Then, the $L^{\prime}_{k}$ matrix after the
perturbation, denoted by $\tilde{L}^{\prime}_{k}$, is written as

$\tilde{L}^{\prime}_{k}=L^{\prime}_{k}+\epsilon\varphi+O(\varepsilon^{2}),$
(20)

where $\varphi\in\mathbb{R}^{k\times n}$ satisfies

$\left\|\varphi\right\|_{F}\leq\left[\frac{1}{4\lambda_{k}^{3/2}}+\frac{1}{2\lambda_{k}^{1/2}(\lambda_{k}-\lambda_{k+1})}\right]\left\|\Phi\right\|_{F}.$
(21)

Replacing an existing landmark by a new one results in a change in a column in
$\Delta$. Theorem 21 provides an upper bound of the amount of perturbation in
$L^{\prime}_{k}$ in response to such a change in $\Delta$. In the following
Theorem 5.5, we exploit Theorem 21 to bound the amount of dislocation of the
embedding coordinate of an arbitrary data point upon a landmark replacement.
The following Lemmas 22, 5.3, and 5.4 aid the proof of Theorem 5.5. Note that
Theorems A.1 and A.2 also assist minor details in the proof of Theorem 5.5.

######  Lemma 5.2.

Let $(X,d)$ be a metric space and $x,y,z\in X$. Assume that
$d(y,z)<\max\\{d(x,y),d(x,z)\\}$. Then, we obtain

$\left|d^{2}(x,y)-d^{2}(x,z)\right|\leq 2d(x,y)d(y,z)+d^{2}(y,z).$ (22)

###### Proof.

The triangle inequality gives $d(x,z)\leq d(x,y)+d(y,z)$. Therefore,

$d^{2}(x,y)-d^{2}(x,z)\geq-2d(x,y)d(y,z)-d^{2}(y,z).$ (23)

Next, we consider another triangle inequality $d(x,y)\leq d(x,z)+d(y,z)$,
which leads to $d(x,y)-d(y,z)\leq d(x,z)$. Because we assumed that
$d(y,z)<\max\\{d(x,y),d(x,z)\\}$, we obtain $d(y,z)\leq d(x,y)$, without the
loss of generality. Therefore, we obtain
$d^{2}(x,y)-2d(x,y)d(y,z)+d^{2}(y,z)\leq d(x,z)^{2}$, i.e.,

$d^{2}(x,y)-d^{2}(x,z)\leq 2d(x,y)d(y,z)-d^{2}(y,z)\leq
2d(x,y)d(y,z)+d^{2}(y,z).$ (24)

Combination of Eqs. (23) and (24) yields Eq. (22). ∎

######  Lemma 5.3.

Let $(X,d)$ be a metric space. Let $\Delta_{1}$ and $\Delta_{2}$ be the
squared distance matrices associated with the sets of landmarks
$\\{x_{1},\ldots,x_{m-1},x_{m}\\}$ and $\\{x_{1},\ldots,x_{m-1},y\\}$,
respectively, where $x_{1},\ldots,x_{m},y\in X$, and $\varepsilon\equiv
d(x_{m},y)$. We further let $L^{\prime(1)}_{k}$ and $L^{\prime(2)}_{k}$ be
$k\times m$ matrices defined by Eq. (13) and associated with
$\\{x_{1},\ldots,x_{m-1},x_{m}\\}$ and $\\{x_{1},\ldots,x_{m-1},y\\}$,
respectively. Furthermore, let $\lambda_{1}\geq\cdots\geq\lambda_{k}$ be the
eigenvalues associated with $\Delta_{1}$ such that $\lambda_{k}>0$ and
$\lambda_{k}>\lambda_{k+1}$. Then, there is a perturbation
$\Delta_{2}=\Delta_{1}+\varepsilon\Phi+O(\varepsilon^{2})$ for small
$\varepsilon$ such that $\left\|\Phi\right\|_{F}\leq 2\sqrt{2K(m-1)}$, where
$K=\max\\{d^{2}(x_{1},x_{m}),\ldots,d^{2}(x_{m-1},x_{m})\\}$. Moreover, if
$z\in X\setminus\\{x_{1},\ldots,x_{m},y\\}$ and $\varepsilon\leq d(z,y)$, then
we obtain

$\left|L^{\prime(2)}_{k}\delta^{(2)}_{z}-L^{\prime(1)}_{k}\delta^{(1)}_{z}\right|_{2}\leq
2\varepsilon\left[d(z,y)\left\|L^{\prime(2)}_{k}\right\|_{2}+\widetilde{K}\sqrt{2K(m-1)}\left|\delta^{(1)}_{z}\right|_{2}\right]+O(\varepsilon^{2}),$
(25)

where $\delta^{(1)}_{z}=\begin{bmatrix}d^{2}(z,x_{1})\\\ \vdots\\\
d^{2}(z,x_{m-1})\\\ d^{2}(z,x_{m})]\end{bmatrix}$,
$\delta^{(2)}_{z}=\begin{bmatrix}d^{2}(z,x_{1})\\\ \vdots\\\
d^{2}(z,x_{m-1})\\\ d^{2}(z,y)\end{bmatrix}$, and
$\widetilde{K}=\frac{1}{4\lambda_{k}^{3/2}}+\frac{1}{2\lambda_{k}^{1/2}(\lambda_{k}-\lambda_{k+1})}$.
Here, $\left\|\cdot\right\|_{2}$ denotes the matrix 2-norm, which we obtain by
substituting
$\left(X_{1},\left\|\cdot\right\|_{X_{1}}\right)=\left(\mathbb{R}^{m},\left|\cdot\right|_{2}\right)$
and
$\left(X_{2},\left\|\cdot\right\|_{X_{2}}\right)=\left(\mathbb{R}^{k},\left|\cdot\right|_{2}\right)$
in Definition 52 in Appendix A.

###### Proof.

First, we obtain

$\Delta_{2}=\Delta_{1}+\left[\begin{array}[]{c|c}\bm{0}_{(m-1)\times(m-1)}&u\\\
\hline\cr u^{\top}&0\end{array}\right],$ (26)

where $\bm{0}_{(m-1)\times(m-1)}$ is the zero matrix of size
$(m-1)\times(m-1)$, and

$u=\begin{bmatrix}d^{2}(x_{1},y)-d^{2}(x_{1},x_{m})\\\
d^{2}(x_{2},y)-d^{2}(x_{2},x_{m})\\\ \vdots\\\
d^{2}(x_{m-1},y)-d^{2}(x_{m-1},x_{m})\end{bmatrix}.$ (27)

For each $r\in\\{1,\ldots,m-1\\}$, we map
$x_{r}\mapsto\tilde{x}_{r},x_{m}\mapsto\tilde{x}_{m}$ and $y\mapsto\tilde{y}$,
where $\tilde{x}_{r},\tilde{x}_{m},\tilde{y}$ are vertices of a triangle,
whose sides are of lengths $d(x_{r},x_{m}),d(x_{r},y)$, and $\varepsilon$, in
$\mathbb{R}^{2}$. This is doable because the metric $d$ satisfies the triangle
inequality. By the law of cosine, we obtain

$d^{2}(x_{r},y)-d^{2}(x_{r},x_{m})=-2d(x_{r},x_{m})\cos(\theta_{r})\varepsilon+\varepsilon^{2},$
(28)

where $\theta_{r}$ is the angle between
$\overrightarrow{\tilde{x}_{r}\tilde{x}_{m}}$ and
$\overrightarrow{\tilde{x}_{m}\tilde{y}}$. By substituting Eq. (28) in Eq.
(27), we obtain

$u=\begin{bmatrix}\phi_{1}\\\ \vdots\\\
\phi_{m-1}\end{bmatrix}\varepsilon+O(\varepsilon^{2}),$ (29)

where

$\phi_{r}\equiv-2d(x_{r},x_{m})\cos(\theta_{r}),\quad r\in\\{1,\ldots,m-1\\}.$
(30)

Then, we define

$\Phi=\begin{bmatrix}0&0&\ldots&0&\phi_{1}\\\ 0&0&\ldots&0&\phi_{2}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\dots&0&\phi_{m-1}\\\
\phi_{1}&\phi_{2}&\ldots&\phi_{m-1}&0\end{bmatrix}.$ (31)

Substituting Eqs. (29), (31) into Eq. (26) yields

$\Delta_{2}=\Delta_{1}+\varepsilon\Phi+O(\varepsilon^{2}).$ (32)

By substituting
$\left(X_{1},\left\langle\cdot,\cdot\right\rangle_{X_{1}}\right)=\left(\mathbb{R}^{m},\left\langle\cdot,\cdot\right\rangle_{E}\right)$
and
$\left(X_{2},\left\langle\cdot,\cdot\right\rangle_{X_{2}}\right)=\left(\mathbb{R}^{k},\left\langle\cdot,\cdot\right\rangle_{E}\right)$
in Theorem A.2, where $\left\langle\cdot,\cdot\right\rangle_{E}$ denotes the
Euclidean inner product, it follows that $L^{\prime(2)}_{k}-L^{\prime(1)}_{k}$
and $L^{\prime(2)}_{k}$ are bounded. Using Theorem A.1, we obtain

$\displaystyle\left|L^{\prime(2)}_{k}\delta^{(2)}_{z}-L^{\prime(1)}_{k}\delta^{(1)}_{z}\right|_{2}$
$\displaystyle=\left|L^{\prime(2)}_{k}\delta^{(2)}_{z}-L^{\prime(2)}_{k}\delta^{(1)}_{z}+L^{\prime(2)}_{k}\delta^{(1)}_{z}-L^{\prime(1)}_{k}\delta^{(1)}_{z}\right|_{2}$
$\displaystyle\leq\left|L^{\prime(2)}_{k}\delta^{(2)}_{z}-L^{\prime(2)}_{k}\delta^{(1)}_{z}\right|_{2}+\left|L^{\prime(2)}_{k}\delta^{(1)}_{z}-L^{\prime(1)}_{k}\delta^{(1)}_{z}\right|_{2}$
$\displaystyle\leq\left\|L^{\prime(2)}_{k}\right\|_{2}\left|\delta^{(2)}_{z}-\delta^{(1)}_{z}\right|_{2}+\left\|L^{\prime(2)}_{k}-L^{\prime(1)}_{k}\right\|_{2}\left|\delta^{(1)}_{z}\right|_{2}$
$\displaystyle\leq\left\|L^{\prime(2)}_{k}\right\|_{2}\left|\delta^{(2)}_{z}-\delta^{(1)}_{z}\right|_{2}+\left\|L^{\prime(2)}_{k}-L^{\prime(1)}_{k}\right\|_{F}\left|\delta^{(1)}_{z}\right|_{2}$
$\displaystyle=\left\|L^{\prime(2)}_{k}\right\|_{2}\left|\begin{bmatrix}0\\\
\vdots\\\ 0\\\
d^{2}(z,y)-d^{2}(z,x_{m})\end{bmatrix}\right|_{2}+\left\|L^{\prime(2)}_{k}-L^{\prime(1)}_{k}\right\|_{F}\left|\delta^{(1)}_{z}\right|_{2}.$
(33)

Lemma 22 implies

$\left|d^{2}(z,y)-d^{2}(z,x_{m})\right|\leq
2d(z,y)\varepsilon+O(\varepsilon^{2}).$ (34)

Because $\Delta_{2}=\Delta_{1}+\varepsilon\Phi+O(\varepsilon^{2})$, by Theorem
21, we can write

$L^{\prime(2)}_{k}=L^{\prime(1)}_{k}+\varepsilon\varphi+O(\varepsilon^{2}),$
(35)

where

$\left\|\varphi\right\|_{F}\leq\widetilde{K}\left\|\Phi\right\|_{F}.$ (36)

By substituting Eqs. (34), (35), and (36) in Eq. (33), we obtain

$\displaystyle(\text{RHS of
Eq.~{}}\eqref{eq:embedding_main_bound_first_part_subs})\leq$
$\displaystyle\left\|L^{\prime(2)}_{k}\right\|_{2}\left[2d(z,y)\varepsilon+O(\varepsilon^{2})\right]+\varepsilon\left\|\varphi\right\|_{F}\left|\delta^{(1)}_{z}\right|_{2}+O(\varepsilon^{2})$
$\displaystyle\leq$
$\displaystyle\varepsilon\left[2d(z,y)\left\|L^{\prime(2)}_{k}\right\|_{2}+\widetilde{K}\left\|\Phi\right\|_{F}\left|\delta^{(1)}_{z}\right|_{2}\right]+O(\varepsilon^{2}).$
(37)

To bound $\left\|\Phi\right\|_{F}$, we use the fact that $\phi^{2}_{r}\leq
4d^{2}(x_{r},x_{m})$, which obeys from Eq. (30), to proceed as follows:

$\left\|\Phi\right\|_{F}=\sqrt{\sum_{r=1}^{m-1}2\phi^{2}_{r}}\leq\sqrt{\sum_{r=1}^{m-1}8d^{2}(x_{r},x_{m})}\leq\sqrt{\sum_{r=1}^{m-1}8K}=2\sqrt{2K(m-1)},$
(38)

where $K$ is defined in the statement of the present Lemma. Using Eq. (38), we
obtain

$\displaystyle\left|L^{\prime(2)}_{k}\delta^{(2)}_{z}-L^{\prime(1)}_{k}\delta^{(1)}_{z}\right|_{2}$
$\displaystyle\leq(\text{RHS of
Eq.~{}}\eqref{eq:embedding_main_bound_first_part_subs})$
$\displaystyle\leq(\text{RHS of
Eq.~{}}\eqref{eq:before_final_bound_first_part})$ $\displaystyle\leq
2\varepsilon\left[d(z,y)\left\|L^{\prime(2)}_{k}\right\|_{2}+\widetilde{K}\sqrt{2K(m-1)}\left|\delta^{(1)}_{z}\right|_{2}\right]+O(\varepsilon^{2}).$
(39)

∎

######  Lemma 5.4.

Let $(X,d)$ be a metric space, and let $\Delta_{1}=[d^{(1)}_{rs}]$ and
$\Delta_{2}=[d^{(2)}_{rs}]$ be the squared distance matrices for
$\\{x_{1},\ldots,x_{m}\\}$ and $\\{x_{1},\ldots,x_{m-1},y\\}$, respectively,
where $x_{1},\ldots,x_{m},y\in X$, and $\varepsilon\equiv d(x_{m},y)$. Then,
we obtain

$\left|\delta^{(2)}_{\mu}-\delta^{(1)}_{\mu}\right|_{2}\leq\frac{4(m-1)\sqrt{K}}{m}\varepsilon+O(\varepsilon^{2})$
(40)

for small $\varepsilon$, where
$\delta^{(1)}_{\mu}=\frac{1}{m}\begin{bmatrix}\sum_{s=1}^{m}d^{(1)}_{1s}\\\
\vdots\\\ \sum_{s=1}^{m}d^{(1)}_{ms}\end{bmatrix}$ and
$\delta^{(2)}_{\mu}=\frac{1}{m}\begin{bmatrix}\sum_{s=1}^{m}d^{(2)}_{1s}\\\
\vdots\\\ \sum_{r=1}^{m}d^{(2)}_{ms}\end{bmatrix}$.

###### Proof.

By combining Eqs. (31) and (32), we obtain

$\Delta_{2}-\Delta_{1}=\varepsilon\begin{bmatrix}0&0&\cdots&0&\phi_{1}\\\
0&0&\cdots&0&\phi_{2}\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\cdots&0&\phi_{m-1}\\\
\phi_{1}&\phi_{2}&\cdots&\phi_{m-1}&0\end{bmatrix}+O(\varepsilon^{2})=\varepsilon\Phi+O(\varepsilon^{2}).$
(41)

Equation (41) implies that

$d^{(2)}_{rs}-d^{(1)}_{rs}=\begin{cases}0,&r,s<m\text{ or }r=s,\\\
\varepsilon\phi_{r}+O(\varepsilon^{2}),&r\in\\{1,\ldots,m-1\\},s=m,\\\
\varepsilon\phi_{s}+O(\varepsilon^{2}),&r=m,s\in\\{1,\ldots,m-1\\}.\end{cases}$
(42)

Therefore,

$\delta^{(2)}_{\mu}-\delta^{(1)}_{\mu}=\frac{1}{m}\begin{bmatrix}\sum_{s=1}^{m}\left(d^{(2)}_{1s}-d^{(1)}_{1s}\right)\\\
\sum_{s=1}^{m}\left(d^{(2)}_{2s}-d^{(1)}_{2s}\right)\\\ \vdots\\\
\sum_{s=1}^{m}\left(d^{(2)}_{m-1,s}-d^{(1)}_{m-1,s}\right)\\\
\sum_{s=1}^{m}\left(d^{(2)}_{ms}-d^{(1)}_{ms}\right)\end{bmatrix}=\frac{\varepsilon}{m}\begin{bmatrix}\phi_{1}\\\
\phi_{2}\\\ \vdots\\\ \phi_{m-1}\\\
\sum_{s=1}^{m-1}\phi_{s}\end{bmatrix}+O(\varepsilon^{2}).$ (43)

Using Eq. (43), we obtain

$\displaystyle\left|\delta^{(2)}_{\mu}-\delta^{(1)}_{\mu}\right|_{2}\leq\left|\delta^{(2)}_{\mu}-\delta^{(1)}_{\mu}\right|_{1}$
$\displaystyle\leq\frac{\varepsilon}{m}\sum_{r=1}^{m-1}\left|\phi_{r}\right|+\frac{\varepsilon}{m}\left|\sum_{s=1}^{m-1}\phi_{s}\right|+\left|O(\varepsilon^{2})\right|$
$\displaystyle\leq\frac{\varepsilon}{m}\sum_{r=1}^{m-1}\left|\phi_{r}\right|+\frac{\varepsilon}{m}\sum_{s=1}^{m-1}\left|\phi_{s}\right|+\left|O(\varepsilon^{2})\right|$
$\displaystyle=2\frac{\varepsilon}{m}\left|\left(\phi_{1},\ldots,\phi_{m-1}\right)\right|_{1}+O(\varepsilon^{2})$
$\displaystyle\leq
2\varepsilon\frac{\sqrt{m-1}}{m}\left|\left(\phi_{1},\ldots,\phi_{m-1}\right)\right|_{2}+O(\varepsilon^{2})$
$\displaystyle=2\varepsilon\frac{\sqrt{m-1}}{m}\sqrt{\frac{\left\|\Phi\right\|^{2}_{F}}{2}}+O(\varepsilon^{2})$
$\displaystyle=\varepsilon\frac{\sqrt{2(m-1)}}{m}\left\|\Phi\right\|_{F}+O(\varepsilon^{2})$
$\displaystyle\leq 4\varepsilon(m-1)\frac{\sqrt{K}}{m}+O(\varepsilon^{2}),$
(44)

where $\left|\cdot\right|_{1}$ denotes the $\ell^{1}$-norm, and we used Eq.
(38) to derive the last inequality. ∎

######  Theorem 5.5.

Let $(X,d)$ be a metric space,
$\psi_{1}=-\frac{1}{2}L^{\prime(1)}_{k}\left(\delta^{(1)}_{z}-\delta^{(1)}_{\mu}\right)$
be the LMDS map associated with the set of landmarks
$\\{x_{1},\ldots,x_{m-1},x_{m}\\}\subseteq X$, and
$\psi_{2}=-\frac{1}{2}L^{\prime(2)}_{k}\left(\delta^{(2)}_{z}-\delta^{(2)}_{\mu}\right)$
be the LMDS map associated with the set of landmarks
$\\{x_{1},\ldots,x_{m-1},y\\}\subseteq X$. Furthermore, let
$\lambda_{1}\geq\cdots\geq\lambda_{k}$ be the eigenvalues of $\Delta_{1}$ such
that $\lambda_{k}>0$ and $\lambda_{k}>\lambda_{k+1}$. With any $z\in
X\setminus\\{x_{1},\ldots,x_{m},y\\}$, we obtain

$\left|\psi_{1}(z)-\psi_{2}(z)\right|_{2}\leq\varepsilon\left\\{\left\|L^{\prime(2)}_{k}\right\|_{2}\left[d(z,y)+\frac{2\sqrt{K}(m-1)}{m}\right]+\widetilde{K}\sqrt{2K(m-1)}\left(\left|\delta^{(1)}_{z}\right|_{2}+\left|\delta^{(1)}_{\mu}\right|_{2}\right)\right\\}+O(\varepsilon^{2})$
(45)

for small $\varepsilon\equiv d(x_{m},y)$.

###### Proof.

Let $z\in X\setminus\\{x_{1},\ldots,x_{m},y\\}$. Then,

$\displaystyle\left|\psi_{1}(z)-\psi_{2}(z)\right|_{2}$
$\displaystyle=\left|-\frac{1}{2}L^{\prime(1)}_{k}\left(\delta^{(1)}_{z}-\delta^{(1)}_{\mu}\right)+\frac{1}{2}L^{\prime(2)}_{k}\left(\delta^{(2)}_{z}-\delta^{(2)}_{\mu}\right)\right|_{2}$
$\displaystyle\leq\frac{1}{2}\left[\left|L^{\prime(2)}_{k}\delta^{(2)}_{z}-L^{\prime(1)}_{k}\delta^{(1)}_{z}\right|_{2}+\left|L^{\prime(1)}_{k}\delta^{(1)}_{\mu}-L^{\prime(2)}_{k}\delta^{(2)}_{\mu}\right|_{2}\right].$
(46)

By Lemma 5.3 (i.e., Eq. (25)), the first term of the RHS of Eq. (46) is
bounded as follows:

$\frac{1}{2}\left|L^{\prime(2)}_{k}\delta^{(2)}_{z}-L^{\prime(1)}_{k}\delta^{(1)}_{z}\right|_{2}\leq\varepsilon\left[d(z,y)\left\|L^{\prime(2)}_{k}\right\|_{2}+\widetilde{K}\sqrt{2K(m-1)}\left|\delta^{(1)}_{z}\right|_{2}\right]+O(\varepsilon^{2}).$
(47)

To bound
$\frac{1}{2}\left|L^{\prime(1)}_{k}\delta^{(1)}_{\mu}-L^{\prime(2)}_{k}\delta^{(2)}_{\mu}\right|_{2}$,
we first obtain

$\displaystyle\left|L^{\prime(1)}_{k}\delta^{(1)}_{\mu}-L^{\prime(2)}_{k}\delta^{(2)}_{\mu}\right|_{2}$
$\displaystyle=\left|L^{\prime(1)}_{k}\delta^{(1)}_{\mu}-L^{\prime(2)}_{k}\delta^{(1)}_{\mu}+L^{\prime(2)}_{k}\delta^{(1)}_{\mu}-L^{\prime(2)}_{k}\delta^{(2)}_{\mu}\right|_{2}$
$\displaystyle\leq\left|L^{\prime(1)}_{k}\delta^{(1)}_{\mu}-L^{\prime(2)}_{k}\delta^{(1)}_{\mu}\right|_{2}+\left|L^{\prime(2)}_{k}\delta^{(1)}_{\mu}-L^{\prime(2)}_{k}\delta^{(2)}_{\mu}\right|_{2}$
$\displaystyle\leq\left\|L^{\prime(1)}_{k}-L^{\prime(2)}_{k}\right\|_{2}\left|\delta^{(1)}_{\mu}\right|_{2}+\left\|L^{\prime(2)}_{k}\right\|_{2}\left|\delta^{(2)}_{\mu}-\delta^{(1)}_{\mu}\right|_{2}$
$\displaystyle\leq\left\|L^{\prime(1)}_{k}-L^{\prime(2)}_{k}\right\|_{F}\left|\delta^{(1)}_{\mu}\right|_{2}+\left\|L^{\prime(2)}_{k}\right\|_{2}\left|\delta^{(2)}_{\mu}-\delta^{(1)}_{\mu}\right|_{2}$
$\displaystyle\leq\left\|L^{\prime(1)}_{k}-L^{\prime(2)}_{k}\right\|_{F}\left|\delta^{(1)}_{\mu}\right|_{2}+\left\|L^{\prime(2)}_{k}\right\|_{2}\cdot
4\varepsilon(m-1)\frac{\sqrt{K}}{m}+O(\varepsilon^{2}).$ (48)

To derive the second and the third inequalities in Eq. (48), we used the fact
that $L^{\prime(1)}_{k}-L^{\prime(2)}_{k}$ and $L^{(2)}_{k}$ are bounded (see
Theorem A.2). To derive the last inequality in Eq. (48), we used Lemma 5.4. By
Theorem 21, we can write
$L^{\prime(2)}_{k}-L^{\prime(1)}_{k}=\varepsilon\varphi+O(\varepsilon^{2})$,
where $\left\|\varphi\right\|_{F}\leq\widetilde{K}\left\|\Phi\right\|_{F}$.
Therefore, we obtain

$\displaystyle\left|L^{\prime(1)}_{k}\delta^{(1)}_{\mu}-L^{\prime(2)}_{k}\delta^{(2)}_{\mu}\right|_{2}\leq(\text{RHS
of Eq.~{}}\eqref{eq:embedding_main_bound_second_part})$
$\displaystyle\leq\varepsilon\widetilde{K}\left\|\Phi\right\|_{F}\left|\delta^{(1)}_{\mu}\right|_{2}+\left\|L^{\prime(2)}_{k}\right\|_{2}\cdot
4\varepsilon(m-1)\frac{\sqrt{K}}{m}+O(\varepsilon^{2})$ $\displaystyle\leq
2\varepsilon\sqrt{K}\left[\widetilde{K}\sqrt{2(m-1)}\left|\delta^{(1)}_{\mu}\right|_{2}+\left\|L^{\prime(2)}_{k}\right\|_{2}\frac{2(m-1)}{m}\right]+O(\varepsilon^{2}),$
(49)

where we used Eq (38) to derive the last inequality in Eq. (49).

By substituting Eqs. (47) and (49) in Eq. (46), we obtain

$\displaystyle\left|\psi_{1}(z)-\psi_{2}(z)\right|_{2}$ $\displaystyle\leq$
$\displaystyle\varepsilon\left[d(z,y)\left\|L^{\prime(2)}_{k}\right\|_{2}+\widetilde{K}\sqrt{2K(m-1)}\left|\delta^{(1)}_{z}\right|_{2}+\frac{1}{2}\cdot
2\varepsilon\sqrt{K}\left[\widetilde{K}\sqrt{2m-2}\left|\delta^{(1)}_{\mu}\right|_{2}+\left\|L^{\prime(2)}_{k}\right\|_{2}\frac{2(m-1)}{m}\right]\right]+O(\varepsilon^{2})$
$\displaystyle=$
$\displaystyle\varepsilon\left\\{\left\|L^{\prime(2)}_{k}\right\|_{2}\left[d(z,y)+\frac{2\sqrt{K}(m-1)}{m}\right]+\widetilde{K}\sqrt{2K(m-1)}\left(\left|\delta^{(1)}_{z}\right|_{2}+\left|\delta^{(1)}_{\mu}\right|_{2}\right)\right\\}+O(\varepsilon^{2}),$
(50)

as desired. ∎

##  6 Numerical results

In this section, we demonstrate the proposed landmark replacement algorithm
with three time series data sets, one of which is synthetic and the other two
are empirical. For each data set, we use the LMDS as the out-of-sampling
embedding method to map the time series data into the two-dimensional
Euclidean space. For each data set, we consider three embedding scenarios. In
the first scenario, referred to as offline embedding with initial landmarks,
we use the first $m$ data points as landmarks and never replace them as newer
data points arrive. In the second scenario, referred to as online embedding,
we update landmarks using Algorithm 2 as each new data point arrives. We
initialize the geometric graph for Algorithm 2 with $\varepsilon=10^{-20}$,
with which the first $m$ nodes are isolated in the initial geometric graph. In
the third scenario, referred to as offline embedding with random landmarks,
after all the points have been revealed, we select $m$ data points uniformly
at random as landmarks and embed the entire data set. We note that the first
and third scenarios are offline out-of-sample embedding methods under
different conditions on the available information.

###  6.1 S-curve

The S-curve is a synthetic data set, constructed by the sklearn package in
python [39]. We generate the S-curve in $\mathbb{R}^{3}$ with $10^{3}$ points,
which we visualize in Fig. 2. Because the S-curve is specified by the
Cartesian coordinate system of the two-dimensional S shape on the $(x,z)$
plane and $0\leq y\leq 2$, we linearly ordered the $10^{3}$ points according
to the $(x,z)$ coordinate so that the point nearest to the right-top corner of
the S shape in the $(x,z)$ plane is the first, and that nearest to the left-
bottom corner of the S shape is the last. The color code in the figure is a
guide to the eyes to indicate the ordering of the points.

The black diamonds in Fig. 2 show the final landmarks in each of the three
embedding scenarios. We set $m=100$. Note that, in the offline embedding with
initial landmarks (shown in panel (a)) and the offline embedding with random
landmarks (shown in panel (c)), the landmarks do not change over time. Figure
2(b) indicates that, with the online embedding, many final landmarks are
devoted to early data points while some newer landmarks are scattered over the
middle and lower parts of the S-curve, covering various parts of the S-curve.
The final $\varepsilon=0.94$ for the online embedding.

We show the results of the LMDS into $\mathbb{R}^{2}$ in Fig. 3 for the three
embedding scenarios. Figure 3(a) shows that the shape of the S-curve is not
preserved in the offline embedding with initial landmarks, presumably because
the set of the first $m$ data points is only good at capturing local geometry
of the S-curve, i.e., the beginning part of the S-curve. In contrast, Figs.
3(b) and 3(c) preserve the structure of the original S-curve well. To
quantitatively assess the quality of embedding, we show in Table 1 the values
of $\sigma_{L}$ and $\sigma$ for each scenario; see the rows with $k=2$, where
we remind that $k$ is the embedding dimension. The normalized stress for the
entire data, $\sigma$, is substantially larger for the offline embedding with
initial landmarks than the online embedding and the offline embedding with
random landmarks. The latter two embedding methods yield similar $\sigma$
values. These results support that our online algorithm performs on par with
the offline embedding with random landmarks, which requires the knowledge of
all the data points beforehand, and better than the offline embedding with
initial landmarks. In terms of $\sigma_{L}$, the offline embedding with
initial landmarks is better than the other two embedding algorithms. This is
because the former has all the $m$ points geometrically close to each other on
the approximately planar part of the S-curve, which renders embedding of the
$m$ landmarks an easy task.

To justify the choice of embedding dimension, $k=2$, we also calculated
$\sigma$ and $\sigma_{L}$ for the three embedding scenarios for $k=1$ and
$k=3$. We show the results in Table 1. We find that $\sigma$ (as well as
$\sigma_{L}$) is much larger than an acceptable range ($\sigma\leq 0.15$ [22])
for any of the three embedding schemes with $k=1$. Therefore, using $k=1$ is
not appropriate. In contrast, all the three embedding scenarios yield
$\sigma,\sigma_{L}<0.01$ when $k=3$. In other words, there is little
distortion by LDMS, which is because the original S-curve is defined in the
three-dimensional Euclidean space. Therefore, we consider that embedding
dimension $k=2$ is an adequate challenge for approximate distance-preserving
embedding of the S-curve data.

Figure 2: Three-dimensional S-curves composed of $10^{3}$ points, shown by the circles. (a) Offline embedding with initial landmarks. (b) Online embedding using Algorithm 2. (c) Offline embedding with random landmarks. The color indicates the artificially defined order in which each data point arrives. Each diamond represents a landmark. In (b), the diamonds are the final landmarks. We set $m=100$. Figure 3: Two-dimensional embedding of the S-curve shown in Fig. 2. (a) Offline embedding with initial landmarks. (b) Online embedding using Algorithm 2. (c) Offline embedding with random landmarks. Each diamond represents a landmark. In (b), the diamonds are the final landmarks. Table 1: Goodness of the LMDS for the three data sets. Columns “Initial” and “Random” represent the offline embedding with the initial and random landmarks, respectively. Data | $m$ | $k$ | Index |  Initial |  Online |  Random
---|---|---|---|---|---|---
S-curve | 100 | 1 | $\sigma_{L}$ |  $0.21$ |  0.40 |  0.29
$\sigma$ |  0.74 |  0.33 |  0.31
2 | $\sigma_{L}$ |  $<0.01$ |  0.15 |  0.13
$\sigma$ |  0.23 |  0.16 |  0.15
3 | $\sigma_{L}$ |  $<0.01$ |  $<0.01$ |  $<0.01$
$\sigma$ |  $<0.01$ |  $<0.01$ |  $<0.01$
SPY | 10 | 1 | $\sigma_{L}$ |  $0.29$ |  0.15 |  0.08
$\sigma$ |  $0.80$ |  0.09 |  0.08
2 | $\sigma_{L}$ |  0.14 |  0.05 |  0.03
$\sigma$ |  0.68 |  0.04 |  0.04
3 | $\sigma_{L}$ |  0.07 |  0.02 |  0.01
$\sigma$ |  0.43 |  0.02 |  0.03
Hospital | 20 | 1 | $\sigma_{L}$ |  0.13 |  0.16 |  0.17
$\sigma$ |  0.61 |  0.22 |  0.23
2 | $\sigma_{L}$ |  $<0.01$ |  0.06 |  0.10
$\sigma$ |  0.40 |  0.11 |  0.11
3 | $\sigma_{L}$ |  $<0.01$ |  $<0.01$ |  0.06
$\sigma$ |  0.29 |  0.07 |  0.07

###  6.2 Stock price time series

Our second data set is the daily adjusted closing prices of the top ten
holding stocks, in terms of market capitalization, from SPDR S&P 500 ETF Trust
(SPY). SPY is a member of Standard & Poor’s depository receipt (SPDR), which
is a family of exchange-traded funds (ETFs). We have normalized the adjusted
closing prices of each stock by the map
$x\mapsto(x-x_{\min})/(x_{\max}-x_{\min})$, where $x$ is the adjusted closing
price of a stock on a valid day, and $x_{\min}$ and $x_{\max}$ are the lowest
and the highest prices of the stock in the observation period, respectively.
As of June 2023, SPY is one of the largest ETFs by market capitalization. The
SPY’s ten largest stock holdings are Apple Inc. (AAPL), Microsoft Corp.
(MSFT), Amazon (AMZN), NVIDIA Corp. (NVDA), Alphabet Inc. Class A (GOOGL),
Tesla Inc. (TSLA), Alphabet Inc. Class C (GOOG), Berkshire Hathaway Inc. Class
B (BRK.B), Meta Platforms Inc. Class A (META), and UnitedHealth Group Inc.
(UNH). We collect the daily adjusted closing prices of these ten stocks from
Yahoo finance from 1/2/2013 to 30/12/2022, resulting in $n=2,518$ valid days
in total.

Each data point, $x_{r}$, where $r\in\\{1,\ldots,n\\}$, is a ten-dimensional
vector, whose each component is each stock’s price on the $r$th valid day. We
use the Euclidean distance to measure the distance between each pair of data
points. Then, we embed $\\{x_{1},\ldots,x_{n}\\}$ into the two-dimensional
Euclidean space. We set the number of landmarks $m=10$. The final
$\varepsilon=0.94$ for the online embedding.

Figure 4: Two-dimensional embedding of the time series of the ten largest
stocks of SPY. Each circle represents the data point on a valid day. The color
of the circle represents the year. Each diamond represents a (final) landmark.
(a) Offline embedding with initial landmarks. (b) Online embedding. (c)
Offline embedding with random landmarks. We set $m=10$. Figure 5: Fraction of
landmarks by year when we embed the time series of the price of the ten
largest stocks of SPY using the online embedding. The time span shown is from
1/15/2013 to 12/30/2022, i.e., from the 11th to the 2,518th days. The dashed
lines represent the boundary between years.

We show in Fig. 4 the embedding coordinates of the daily adjusted closing
prices of ten stocks. Unlike in Fig. 2(b), which shows that a majority of the
final landmarks are data points arriving in early time steps (shown by the
diamonds), seven out of the ten final landmarks shown in Fig. 4(b) are data
points arriving in 2018 or later, and the remaining three landmarks are in
2013. To examine evolution of landmarks over years, we show in Fig. 5 the
fraction of the year of the landmarks at each day between the 11th
($=(m+1)$th) and the last days. Specifically, to calculate the fraction of the
year of the landmarks, we count the number of landmarks that correspond to
each year and divide the number by $m$ $(=10)$. For example, if the fraction
of year 2013 is 0.2 at time $t$, it means that two landmarks at time $t$ are
dated year 2013. Figure 5 indicates that our online embedding algorithm
selects fewer data points in 2016, 2018, and 2019 as landmarks, even
transiently. An intuitive explanation of this phenomenon is that the data
points in these years do not explore new part of the embedding space (see Fig.
4(b)), or presumably in the original ten-dimensional space, either.

The trajectory in the embedding space obtained with the offline embedding with
initial landmarks looks qualitatively different from those obtained with the
online embedding and the offline embedding with random landmarks.
Specifically, the data points in 2022 heavily overlap with those between 2016
and 2021 in the case of the offline embedding with initial landmarks (see Fig.
4(a)). In contrast, with the online embedding and the offline embedding with
random landmarks, the data points in 2022 stretch into part of the embedding
space that has been unexplored by the trajectory in the prior years (see Figs.
4(b) and 4(c)). Similarly to the case of the S-curve data set, we hypothesize
that the set of the first $m$ data points, which are the landmarks in the
offline embedding with initial landmarks, only captures local geometry of the
data set near the first $m$ data points and fails to capture the structure of
the later data points that are far from the first $m$ landmarks. To discuss
this possibility, we show $\sigma$ and $\sigma_{L}$ for the three embedding
scenarios in Table 1. The table indicates that the online embedding and the
offline embedding with random landmarks perform better than the offline
embedding with initial landmarks in terms of both $\sigma$ and $\sigma_{L}$
for this data set, affirming our hypothesis. The table also indicates that the
online embedding and the offline embedding with random landmarks produce
almost the same quality of embedding (i.e., $\sigma=0.04$).

We show in Table 1 the values of $\sigma$ and $\sigma_{L}$ for embedding
dimensions $k=1,3$ and for each of the three embedding scenarios. The offline
embedding with initial landmarks gives large $\sigma$ values, implying poor
embedding, for both $k=1$ and $k=3$, as expected. It gives an acceptable
$\sigma_{L}$ value for $k=3$, which only suggests successful embedding of the
landmarks, not the non-landmark points. With the online embedding and offline
embedding with random landmarks, $\sigma$ is sufficiently small (i.e., $\leq
0.1$) for both $k=1$ and $k=3$ as well as for $k=2$. However, the gap in terms
of the $\sigma$ value is larger between $k=1$ and $k=2$ than between $k=2$ and
$k=3$, and the same holds true for $\sigma_{L}$. We chose $k=2$ for these
reasons.

###  6.3 Hospital

As a third data set, we analyze a temporal network data set using the LMDS-
based embedding method proposed in our previous study [77]. The data set,
which we refer to as the Hospital data set, is a sequence of time-stamped
contact events between two individuals in a hospital ward in Lyon, France,
recorded by the SocioPatterns project [41]. We transform the data set into a
continuous-time temporal networks, i.e., weighted adjacency matrix as a
function of continuous time, using the tie-decay network model [69]. The
constructed temporal network contains 75 nodes, of which 11 nodes represent
medical doctors, 35 nurses, and 29 patients. There are $1,139$ edges and
$32,424$ contacts in total. The data set was recorded between 1:00 PM,
12/6/2010, Monday and 2:00 PM, 12/10/2010, Friday. There are $n=9,453$ time
points at which there is at least one contact event, which we use as $t_{1}$,
$\ldots$, $t_{n}$.

We review tie-decay networks and the embedding method based on them [77] in
Appendix A. The tie-decay network is a $75\times 75$ matrix that varies over
time, and each of its entry represents the time-dependent strength of the tie
between a pair of nodes. The time being continuous implies that the network is
not limited to the times at which at least one contact event occurs. We embed
the tie-decay network obtained from the Hospital data into the two-dimensional
Euclidean space. The embedding method is based on the LMDS and therefore
requires a distance measure between pairs of data points, which are networks
in the present case. We employ the unnormalized Laplacian distance as network
distance measure because it yielded satisfactory embedding in our previous
study [77]. We set $m=20$ and the decay rate of the tie strength
$\alpha=10^{-2}$ (see Appendix A for the definition of $\alpha$). The final
$\varepsilon=22.79$ for the online embedding.

We show in Fig. 6 the hourly averages of the embedding coordinates for each of
the three embedding scenarios. For example, the coordinate labeled 2 PM
represents the averaged coordinate of the trajectory between 2:00 PM and 3:00
PM. Similar to the case of the S-curve and SPY data sets, the offline
embedding with initial landmarks produces an embedding trajectory that is
qualitatively different from the online embedding and the offline embedding
with random landmarks. Specifically, the embedded data points are roughly
linearly distributed in the two-dimensional embedding space (see Fig. 6(a)).
In contrast, the trajectory is U-shaped in the case of the online embedding
(see Fig. 6(b)). The shape of the trajectory is roughly intermediate between
the trajectories for these two cases for the offline embedding with random
landmarks (see Fig. 6(c)).

To illustrate how the landmarks are distributed over time and days, we show in
Fig. 7 the time distribution of the $m=20$ landmarks for each of the three
embedding scenarios. Each diamond represents a landmark. By definition, the
landmarks for the offline embedding with initial landmarks are concentrated
around 1:00 PM, Monday, specifically, between 1:00 PM to 1:12 PM (see Fig.
7(a)). The landmarks for the offline embedding with random landmarks are
widely distributed over the different days and time (Fig. 7(c)), which is also
expected. In contrast, Fig. 7(b) indicates that all but two final landmarks
for the online embedding are the networks between 1:00 PM and 2:00 PM, Monday,
and the other two landmarks are from Monday at 2:14 PM and Tuesday at 10:16
AM. We interpret that this result occurs because the embedding trajectory does
not much explore new part of the embedding space after Monday such that the
landmarks from Monday are mostly sufficient.

Table 1 indicates that the offline embedding with initial landmark yields
$\sigma=0.40$ (see the rows with $k=2$), which is substantially larger than
the $\sigma$ values for the other two embedding scenarios. Therefore, the
offline embedding with initial landmarks is not satisfactory despite that it
is good at capturing the local geometry of the set of the initial landmarks
(with $\sigma_{L}<0.01$). We obtain $\sigma=0.11$ for both of the online
embedding and the offline embedding with random landmarks. These $\sigma$
values are roughly as small as that when we use all the 9,453 data points at
which any contact event arrives as landmarks (i.e., $\sigma=0.09$ [77]).
Therefore, we conclude that our online embedding algorithm with a small number
of landmarks (i.e., $m=20$) and the two-dimensional embedding space is fairly
satisfactory for this data set.

We previously showed for this data set that two networks at the same time of
different days tend to be more similar than two networks whose times of the
day are different [77]. To examine this property for our online landmark
replacement algorithm, we compute the distances between arbitrary pairs of
hours in each of Figs. 6(a), 6(b), and 6(c). For each embedding scenario, we
group the distance values into two sets. The first set is composed of the
distance values for an arbitrary pair of networks at the same time of the day.
The second set is composed of the distance values for an arbitrary pair of
networks at different times of the day. If the trajectory is similar across
the different days, then the distance for the first set will be smaller than
that for the second set. Therefore, we run the Student’s t-test for unequal
sample sizes to test the difference in the average distance between the two
sets of distance values. The $p$-values for the offline embedding with initial
landmarks, online embedding, and offline embedding with random landmarks are
equal to $2.46\times 10^{-9},8.08\times 10^{-15}$, and $5.74\times 10^{-15}$,
respectively. Therefore, the Hospital temporal contact network produces
similar trajectories over the different days, including the time, according to
the t-test, no matter which of the three embedding scenarios we adopt.
However, the $p$-value for the offline embedding with initial landmarks is
roughly $10^{6}$ times larger than those for the online embedding and the
offline embedding with random landmarks. Furthermore, the two-dimensional
embedding in which we use all the $n=9,453$ time points as landmarks yields
$p=3.67\times 10^{-15}$, which is roughly $10^{6}$ times smaller than that for
the offline embedding with initial landmarks, but only roughly $2.2$ and $1.6$
times smaller than those for the online embedding and the offline embedding
with random landmarks, respectively. This result lends another support to the
strength of our online embedding algorithm.

Last, we highlight that Table 1 also provides the $\sigma$ and $\sigma_{L}$
values for embedding dimensions $k=1$ and $k=3$. Regardless of the $k$
($\in\\{1,2,3\\}$) value, $\sigma$ for the offline embedding with initial
landmarks is unacceptably large and larger than $\sigma$ for the online
embedding and the offline embedding with random landmarks. For the online
embedding and the offline embedding with random landmarks with $k=1$, we
obtain $\sigma=0.22$ and $\sigma=0.23$, respectively. These $\sigma$ values
are higher than a guideline cut-off value of $0.15$ [22]. In contrast,
$\sigma$ is small enough (i.e., $\sigma\leq 0.11$) for these two embedding
algorithms when $k=2$ and $k=3$. The gap in $\sigma$ is larger between $k=1$
and $k=2$ than between $k=2$ and $k=3$, which is why we have mainly used
$k=2$. These results on the dependence of $\sigma$ on $k$ are similar to those
for the SPY data.

Figure 6: Hourly average embedding coordinate of the Hospital temporal network
across weekdays. (a) Offline embedding with initial landmarks. (b) Online
embedding. (c) Offline embedding with random landmarks. We set $m=20$ and
$\alpha=10^{-2}$. Figure 7: Distribution of the final landmarks over time for
the Hospital temporal network. The midnight corresponds to $0$ and $24$ on the
horizontal axis. (a) Offline embedding with initial landmarks. (b) Online
embedding. (c) Offline embedding with random landmarks.

###  6.4 Computation time

In this section, we numerically examine the computation time and time courses
of $\varepsilon$ for our online embedding algorithm. We set the embedding
dimension $k=2$. We carried out the computations on a Macbook pro with an
Apple M1 chip as the processing unit and 8 GB of RAM.

In practice, a landmark replacement does not occur every time a new data
point, $x_{r}$, arrives. Therefore, the actual computation time may be smaller
than the theoretical results derived in Corollary 3.8, i.e., $O(n^{3})$, which
corresponds to the worst case. We show in Figs. 8(a), 8(b), and 8(c) the time
consumed for sequentially embedding the first $n$ data points in the S-curve,
SPY, and Hospital data sets, respectively. Note that we used $n$ to denote the
number of all data points in the previous sections, whereas here we use $n$ to
denote an arbitrary number of data points. However, we believe that there is
no confusion.

Figures 8(a) and 8(b) suggest that the actual computation time approximately
scales as $O(n^{2})$ rather than $O(n^{3})$. In the figure, we show $\propto
n^{2}$ and $\propto n^{3}$ lines as a guide to the eyes, where $\propto$
represents “in propotion to”. We need to be careful when assessing the
computation time for the Hospital data because it is a temporal network,
whereas the S-curve and SPY data are conventional multidimensional time
series. We recall that we have used the unnormalized Laplacian network
distance measure for the Hospital data and the Euclidean distance for the
S-curve and SPY data. Computing the former takes much longer time,
specifically, $O(n^{3})$ [77], than computing the latter because the former
requires computation of eigenvalues. Therefore, we precompute the Laplacian
distance for each pair of static networks in the tie-decay temporal network.
Then, we carry out our online embedding algorithm by referring to the
precomputed (and stored) Laplacian distance values. We show in Fig. 8(c) two
types of computation time for embedding the Hospital data. One is the total
computation time. The other is the computation time excluding the calculation
of the Laplacian distance between each pair of static networks. We find that
the computation time without the calculation of the Laplacian distance roughly
scales as $O(n^{3})$. Therefore, the calculation of the Laplacian distance is
not a bottleneck of computation. The total computation time starts at a
relatively large value but seems to approach the $O(n^{3})$ scaling as $n$
increases. In conclusion, the numerically evaluated time complexity is at most
cubic and smaller for the S-curve and SPY data. These results are consistent
with our theoretical results.

Figure 8: Numerically evaluated computation time of the proposed online
embedding algorithm. (a) S-curve. (b) SPY. (c) Hospital. The green and orange
straight lines represent $\propto n^{2}$ and $\propto n^{3}$, respectively.
The intercept of these lines is arbitrary. In (c), we show the total
computation time and the computation time excluding the computation of the
Laplacian distance matrix by different lines.

##  7 Discussion

We developed an algorithm of replacing landmarks in an online fashion for out-
of-sample embedding techniques. In particular, we employed the framework of
geometric graphs and dominating sets to capture geometry of landmarks and
other data points. We then provided mathematical underpinning of the proposed
algorithm on the computation time and perturbation of the embedding coordinate
upon a landmark replacement event. We finally carried out numerical
simulations to verify the performance of the algorithm combined with the LMDS.
On the three data sets examined, an embedding dimension of two was practically
enough, and the algorithm performed on par with the offline counterpart of the
LMDS in which the fixed landmarks are selected uniformly at random from all
the data points (which is by definition impossible for any online algorithm).

Our online landmark replacement scheme crucially depends on variable
$\varepsilon$, including its initial value. Variable $\varepsilon$ regulates
the number of edges in the geometric graph used in Algorithm 2. If the arrival
of a new data point triggers an actual landmark replacement event, then
$\varepsilon$ increases. We numerically verified that $\varepsilon$
drastically increases when the trajectory explores new part of the data space
and that $\varepsilon$ saturates or only gradually increases after some time
(see Appendix A). Our algorithm forces $\varepsilon$ to grow monotonically in
time and lacks a mechanism to limit and shrink the size of $\varepsilon$. A
large $\varepsilon$ implies that the update of the geometric graph is
insensitive to the position of new data points. For example, if $\varepsilon$
is extremely large, then any new data point would be connected to all the
existing nodes. We suspect that $\varepsilon$ in later stages of our
simulations is larger than would be necessary such that the minimal dominating
set of the geometric graph tends to contain less than $m$ nodes and therefore
our landmark choice is suboptimal. If $\varepsilon$ were smaller such that the
minimal dominating set contains just $m$ nodes, then these $m$ nodes may be a
better set of $m$ landmarks. To realize such an online embedding scheme
requires a mechanism to both increase and decrease $\varepsilon$ without a
formidable computational cost. Extension of our algorithm to this direction is
left for future work.

We theoretically showed that the time complexity of the proposed algorithm is
$O(n^{3})$. Numerical experiments suggested $O(n^{2})$ or $O(n^{3})$ time
depending on the data set. If we run the LMDS every time a new data point
$x_{m+r}$ arrives, with $\\{x_{1},\ldots,x_{m+r}\\}$ as landmarks, to generate
an embedding trajectory [77], then the total time complexity for adding the
$(m+1)$th to the $(m+n)$th data points is
$\sum_{r=1}^{n}O\left((m+r)^{3}\right)=O(n^{4})$. Therefore, our online
algorithm reduces the computational time by a factor of $n$. Nevertheless,
faster online algorithms are desirable. Because our landmark replacement
algorithm just uses the distance matrix as input, it does not have to be
combined with LMDS. Among the online and out-of-sample algorithms discussed in
section 1, L-ISOMAP, the online L-ISOMAP, and the landmark diffusion maps can
be combined with our proposed landmark replacement algorithm. This is because
these three algorithms use a fixed number of landmarks. While L-ISOMAP shares
the same time complexity as LMDS [20], the online L-ISOMAP [25] and the
landmark diffusion maps [63] are faster than $O(n^{3})$. We used LMDS only
because Theorem 21, a theorem associated with LMDS, allows us derive
mathematical properties (Lemma 5.3 and Theorem 5.5) of the embedding using
LMDS. Exploring online landmark replacement with the online L-ISOMAP, landmark
diffusion maps, and similar algorithms may be interesting.

A potential for future work is to incorporate change-point detection into the
landmark replacement framework to suppress the computation time. In the
change-point detection problem, we examine the presence of a change in a given
time series at each time step, usually via a statistical test [32, 45, 67, 26,
27, 61]. One idea is to trigger the landmark replacement algorithm only when a
new data point is detected to mark a change point. In this manner, we can
reduce the number of times that one calls the actual landmark replacement
routine. A potential challenge is how to maintain $L$ to be the dominating set
when a new data point is an outlier to an intermediate extent such that it is
not detected to be a change point whereas it is an isolated node in the
geometric graph when added. Devising landmark replacement algorithms with
change-point detection is worth exploring.

##  8 Acknowledgments

We thank Lorenzo Livi for discussion. The work of Naoki Masuda was supported
in part by the Air Force Office of Scientific Research (AFOSR) European Office
under Grant FA9550-19-1-7024, in part by the National Science Foundation or
NSF under grant DMS-2052720, in part by JSPS KAKENHI under grants JP21H04595
and 23H03414, and in part by the Japan Science and Technology Agency (JST)
under grant JPMJMS2021.

##  Appendix A Linear operators

In this section, we review basic linear operator theory that we have used in
the main text.

######  Definition 3.

Consider normed spaces $\left(X_{1},\left\|\cdot\right\|_{X_{1}}\right)$ and
$\left(X_{2},\left\|\cdot\right\|_{X_{2}}\right)$, and let $T:X_{1}\rightarrow
X_{2}$ be a linear map. Then, $T$ is bounded if there exists
$\widetilde{M}\geq 0$ such that

$\left\|T(x)\right\|_{X_{2}}\leq\widetilde{M}\left\|x\right\|_{X_{1}}$ (51)

for all $x\in X_{1}$.

######  Definition 4.

Consider normed spaces $\left(X_{1},\left\|\cdot\right\|_{X_{1}}\right)$ and
$\left(X_{2},\left\|\cdot\right\|_{X_{2}}\right)$, and let $T:X_{1}\rightarrow
X_{2}$ be a bounded linear map. We define the operator norm of $T$ as follows:

$\left\|T\right\|_{\textnormal{op}}=\inf\\{\widetilde{M}\geq
0:\left\|T(x)\right\|_{X_{2}}\leq\widetilde{M}\left\|x\right\|_{X_{1}}\forall
x\in X_{1}\\}.$ (52)

Note that Definitions 3 and 52 depend on the choice of norms.

######  Theorem A.1.

(Einsiedler and Wart [57]) Let
$\left(X_{1},\left\|\cdot\right\|_{X_{1}}\right)$ and
$\left(X_{2},\left\|\cdot\right\|_{X_{2}}\right)$ be normed spaces, and
$T:X_{1}\rightarrow X_{2}$ be a bounded linear map. Then,

$\left\|T(x)\right\|_{X_{2}}\leq\left\|T\right\|_{\textnormal{op}}\left\|x\right\|_{X_{1}},$
(53)

for all $x\in X_{1}$.

Recall that a metric space is complete if every Cauchy sequence converges to
an element in the metric space, and that an inner product space
$\left(X,\left\langle\cdot,\cdot\right\rangle\right)$ is a Hilbert space if it
is complete with respect to its induced norm
$\left\|\cdot\right\|=\sqrt{\left\langle\cdot,\cdot\right\rangle}$.

######  Theorem A.2.

(Halmos [38]) Let
$\left(X_{1},\left\langle\cdot,\cdot\right\rangle_{X_{1}}\right)$ and
$\left(X_{2},\left\langle\cdot,\cdot\right\rangle_{X_{2}}\right)$ be Hilbert
spaces. If $X_{1}$ is finite-dimensional and $T:X_{1}\rightarrow X_{2}$ is
linear, then $T$ is bounded.

##  Appendix A Tie-decay network and its embedding

In this section, we briefly review tie-decay networks and its embedding.

Denote the number of nodes of the temporal network by $N$. Consider a sequence
of time-stamped contact events each of which occurs between a pair of nodes.
One can represent each contact event as a tuple $(i_{e},j_{e},t_{e})$, where
$i_{e}$ and $j_{e}$ are nodes, and $t_{e}$ $(\geq 0)$ is the time of the
contact event. For the Hospital data set, each contact event is undirected
such that the order of $i_{e}$ and $j_{e}$ does not matter. Let
$B(t)=[b_{ij}(t)]$ be a time-dependent $N\times N$ weighted adjacency matrix
of the tie-decay network, where $b_{ij}(t)$ represents the weight of edge
$(i,j)$ at time $t\in\mathbb{R}_{\geq 0}$. For each $(i,j)$, we set

$b_{ij}(t)=\sum_{r;\tilde{t}_{r}\leq
t}e^{-\alpha(t-\tilde{t}_{r})}H(t-\tilde{t}_{r}),$ (54)

where $H$ is the piecewise Heaviside step function defined by $H(x)=1$ for
$x\geq 0$, and $H(x)=0$ for $x<0$, and $\tilde{t}_{r}$ is the time of the
$r$th event on edge $(i,j)$, and $\alpha$ is a parameter [69, 73, 68, 75]. An
event on $(i,j)$ increases $b_{ij}(t)$ by $1$. We set $\alpha=10^{-2}$ in our
numerical simulations.

The exponentially decaying nature of Eq. (54) is convenient for updating
matrix $B$ upon any event arrival. To explain this, we consider a sequence of
times $0\leq t_{1}<t_{2}<\cdots<t_{n}$, where $t_{r}$ is the $r$th time at
which at least one event occurs in the entire network. We denote by $A_{r}$
the adjacency matrix of the network composed of all the events at time
$t_{r}$. The input data is then equivalent to the time series of adjacency
matrices $\\{A_{1},A_{2},\ldots\\}$. Then, the weighted adjacency matrix of
the tie-decay network is given by

$B(t)=\sum_{r;t_{r}\leq t}e^{-\alpha(t-t_{r})}A_{r},$ (55)

which allows a convenient updating rule:

$B(t_{r+1})=e^{-\alpha(t_{r+1}-t_{r})}B(t_{r})+A_{r+1}.$ (56)

We embed $\\{B(t_{1}),\ldots,B(t_{n})\\}$ using the LMDS. This method is
advantageous in that one can produce a continuous-time trajectory of a given
tie-decay network, i.e., including the times at which no contact event occurs,
if we use an appropriate distance measure for networks [77]. In the
demonstration with the Hospital data set, we do not exploit this continuous-
time nature of the trajectory because our aim is to assess the performance of
our online landmark replace algorithm. However, we use this method.

To run the algorithm for a given $\alpha$ value, we first construct the tie-
decay network from the input adjacency matrices, $\\{A_{1},\ldots,A_{n}\\}$,
using Eq. (55). Second, we specify the landmarks. In the original method [77],
we use $B(t_{1})$, $\ldots$, $B(t_{n})$ as landmarks and embed $B(t)$, where
$t\notin\\{t_{1},\ldots,t_{n}\\}$, using the LMDS. Here we only use the $m$
landmarks that are a subset of $\\{B(t_{1}),\ldots,B(t_{n})\\}$. Third, we run
the LMDS to embed all $B(t)$’s that are not landmarks. For $B(t_{r})$ that are
not landmarks, the usual procedures of the LMDS determine their embedding
coordinates. To embed $\\{B(t):t\geq 0\\}$, where
$t\notin\\{t_{1},\ldots,t_{n}\\}$, we use the exponentially decaying nature of
the tie-decay network model to obtain

$B(t)=e^{-\alpha(t-t_{r})}B(t_{r}),$ (57)

where $r$ is the unique value of $r\in\\{1,\ldots,n\\}$ that satisfies
$t_{r}<t<t_{r+1}$. Using Eq. (12), we obtain the embedding coordinate for
$B(t)$ for any $t_{r}<t<t_{r+1}$ by

$\psi(B(t))=-\frac{1}{2}L^{\prime}_{k}(\delta_{B}(t)-\delta_{\mu}),$ (58)

where

$\delta_{B}(t)=\begin{bmatrix}d^{2}(B(t),B(t_{1}))\\\ \vdots\\\
d^{2}(B(t),B(t_{n}))\end{bmatrix}.$ (59)

We use the unnormalized Laplacian network distance as $d$. The Laplacian
matrix $\tilde{L}(t)$ of a tie-decay network at time $t$ is given by
$\tilde{L}(t)\equiv\tilde{D}(t)-B(t)$, where $\tilde{D}(t)$ is the diagonal
matrix of which the $i$th diagonal entry is
$\displaystyle\sum_{j=1}^{N}b_{ij}(t)$. The unnormalized Laplacian network
distance is given by

$d_{\tilde{L}}(\tilde{L}(t_{r}),\tilde{L}(t_{s}))=\sqrt{\sum_{k=1}^{N}\left[\rho_{k}(\tilde{L}(t_{r}))-\rho_{k}(\tilde{L}(t_{s}))\right]^{2}}\,,$
(60)

where $\rho_{k}(P)$ is the $k$th smallest eigenvalue of symmetric matrix $P$
[64, 60, 31]. One can replace $\sum_{k=1}^{N}$ by $\sum_{k=2}^{N}$ in Eq. (60)
because the smallest eigenvalue of a Laplacian matrix is always $0$. Using
Eqs. (57), we can simplify (60) as follows:

$d_{\tilde{L}}(\tilde{L}(t),\tilde{L}(t_{s}))=\sqrt{\sum_{k=2}^{N}\left[e^{-\alpha(t-t_{r})}\rho_{k}(\tilde{L}(t_{r}))-\rho_{k}(\tilde{L}(t_{s}))\right]^{2}}\,.$
(61)

##  Appendix A Growth of $\varepsilon$ over time

We plot in Figs. 9(a), 9(b), and 9(c) the values of $\varepsilon$ against the
number of data points for the S-curve, SPY, and Hospital data, respectively.
In Fig. 9(a), $\varepsilon$ rapidly grows till approximately the first 200
data points. In Fig. 9(b), $\varepsilon$ starts growing faster than before at
$r\approx 1,500$ because of the nature of the SPY data; the stock prices
started to explore different part of the space of the data, roughly in the
beginning of year 2020 (i.e., 1,761th day). Although this interpretation is
consistent with Fig. 4(b), we further support this claim by the following
computation. For simplicity, we calculated the average of the ten-dimensional
vector whose each entry is the price of a stock averaged over all valid days
of a year. The distance between the thus obtained vectors of two years, which
we averaged over all pairs of years between 2013 and 2019, was equal to
$1.80\times 10^{-4}$. In contrast, the distance between the vector for year
2020 and that for any year between 2013 and 2019 was equal to $4.90\times
10^{-4}$ on average. Therefore, the coordinates between 2013 and 2019 in the
original space stay closer to each other than the coordinates in 2020. In Fig.
9(c), the value of $\varepsilon$ suddenly increases by a large amount on
Monday at 2:08 PM, Monday at 2:14 PM, and Tuesday at 10:16 AM. Note that we
use the logarithmic scale of $n$ in Fig. 9(c) to uncover the first two big
leaps in $\varepsilon$ values.

Figure 9: Evolution of $\varepsilon$ as a function of time. (a) S-curve. (b)
SPY. (c) Hospital.

This last phenomenon is considered to be because the embedding trajectory from
Monday between 2:00 PM and 3:00 PM and Tuesday between 10:00 AM and 11:00 AM
explores unexplored part of the data space compared to other segments of the
same trajectory in the past. As evidence, we show in Fig. 10 fifteen one-hour
segments of the trajectory of the Hospital temporal network, one per panel.
The segments of the trajectory shown in Fig. 10 that include Monday 2:08 AM
and 2:14 PM explore new part of the space relative to the segments in the past
(i.e., the segment between Monday 1:00 PM and 2:00 PM). Similarly, apart from
the segment containing Tuesday 10:16 AM, all segments up to Tuesday 10:00 AM
are confined below $y=6$ and to the right to $x=-25$. After the third sudden
increment, $\varepsilon$ $(=22.79)$ is large enough that the existing
landmarks cover the incoming networks. We recall that the $x$ and $y$ ranges
shown in Fig. 6(b) are appoximately $x\in(-15,10)$ and $y\in(-2.5,5.5)$,
respectively, which are a small subset of the embedding space shown in each
panel of Fig. 10. This is because we show the hourly averages of the embedding
coordinate in Fig. 6, while we show the embedding coordinate at individual
times in Fig. 10; the latter is not an averaged quantity and therefore
accompanies larger fluctuations.

Figure 10: Two-dimensional trajectory of the Hospital temporal network in each
hour from Monday 1:00 PM to 8:00 PM and from Tuesday 6:00 AM to 3:00 PM.
Because 95.6% of the contacts occur between 6:00 AM to 8:00 PM throughout the
weekdays, we have omitted the trajectories outside these hours. We show the
trajectory on Monday and Tuesday in different colors.

While there is no control or a stopping criterion for $\varepsilon$, it should
stop growing after its value has become sufficiently large because any newly
arriving data points, including outliers, would be connected to an existing
landmark when $\varepsilon$ is sufficiently large. In fact, Fig. 9(c) shows
that $\varepsilon$ saturates after Tuesday at 10:16 AM (i.e., after $n=1,838$
in the figure). Although $\varepsilon$ does not saturate for the other two
data sets (see Figs. 9(a) and 9(b)), we expect that it roughly saturates if
the data are longer in the case of the SPY data. In the case of the S-curve,
$\varepsilon$ would continue to increase if the S-like shape further
continues. However, this situation is obviously artificial.

## References

  * Kruskal, [1964] Kruskal, J. B.  Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis.  _Psychometrika_ , 29:1–27, 1964.
  * Hedetniemi and Laskar, [1990] Hedetniemi, S. T. and Laskar, R. C.  Bibliography on domination in graphs and some basic definitions of domination parameters.  _Discrete Mathematics_ , 86:257–277, 1990.
  * Parekh, [1991] Parekh, A. K.  Analysis of a greedy heuristic for finding small dominating sets in graphs.  _Information Processing Letters_ , 39:237–240, 1991.
  * Carreira-Perpiñán, [1997] Carreira-Perpiñán, M. Á.  A review of dimension reduction techniques.  Technical report, Department of Computer Science, University of Sheffield, 1997.
  * King and Tzeng, [1997] King, G.-H. and Tzeng, W.-G.  On-line algorithms for the dominating set problem.  _Information Processing Letters_ , 61:11–14, 1997.
  * Graepel et al., [1998] Graepel, T., Herbrich, R., Bollmann-Sdorra, P., and Obermayer, K.  Classification on pairwise proximity data.  In _Proc. Advances in Neural Information Processing Systems_ , pages 438–444. 1998.
  * Hall et al., [1998] Hall, P. M., Marshall, A. D., and Martin, R. R.  Incremental eigenanalysis for classification.  In _Proc. the British Machine Vision Conference_, pages 29.1–29.10. 1998.
  * Knuth, [1998] Knuth, D. E.  _The Art of Computer Programming: Volume 3: Sorting and Searching_.  Addison-Wesley, Boston, MA, 2nd edition, 1998.
  * Dasgupta, [2000] Dasgupta, S.  Experiments with random projection.  In _Proc. the 16th Conference on Uncertainty in Artificial Intelligence_ , pages 143–151. 2000.
  * Roweis and Saul, [2000] Roweis, S. T. and Saul, L. K.  Nonlinear dimensionality reduction by locally linear embedding.  _Science_ , 290:2323–2326, 2000.
  * Tenenbaum et al., [2000] Tenenbaum, J. B., De Silva, V., and Langford, J. C.  A global geometric framework for nonlinear dimensionality reduction.  _Science_ , 290:2319–2323, 2000.
  * Bingham and Mannila, [2001] Bingham, E. and Mannila, H.  Random projection in dimensionality reduction: applications to image and text data.  In _Proc. the 7th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 245–250. 2001.
  * Pekalska et al., [2001] Pekalska, E., Paclik, P., and Duin, R. P.  A generalized kernel approach to dissimilarity-based classification.  _Journal of Machine Learning Research_ , 2:175–211, 2001.
  * Artac et al., [2002] Artac, M., Jogan, M., and Leonardis, A.  Incremental PCA for on-line visual learning and recognition.  In _Proc. 2002 International Conference on Pattern Recognition_ , volume 3, pages 781–784. 2002.
  * De Silva and Tenenbaum, [2002] De Silva, V. and Tenenbaum, J. B.  Global versus local methods in nonlinear dimensionality reduction.  In _Proc. Advances in Neural Information Processing Systems_ , pages 721–728. 2002.
  * Eidenbenz, [2002] Eidenbenz, S. J.  Online dominating set and variations on restricted graph classes.  Technical report, Department of Computer Science, ETH Zurich, 2002.
  * Achlioptas, [2003] Achlioptas, D.  Database-friendly random projections: Johnson-Lindenstrauss with binary coins.  _Journal of Computer and System Sciences_ , 66:671–687, 2003.
  * Penrose, [2003] Penrose, M.  _Random Geometric Graphs_.  Oxford University Press, New York, NY, 2003.
  * Bengio et al., [2004] Bengio, Y., Paiement, J. F., Vincent, P., Delalleau, O., Le Roux, N., and Ouimet, M.  Out-of-sample extensions for LLE, ISOMAP, MDS, Eigenmaps, and Spectral clustering.  In _Proc. Advances in Neural Information Processing Systems_ , volume 16, pages 177–184. 2004.
  * De Silva and Tenenbaum, [2004] De Silva, V. and Tenenbaum, J. B.  Sparse multidimensional scaling using landmark points.  Technical report, Stanford University, 2004.
  * Law et al., [2004] Law, M. H., Zhang, N., and Jain, A. K.  Nonlinear manifold learning for data stream.  In _Proc. the 2004 SIAM International Conference on Data Mining_ , pages 33–44. 2004.
  * Borg and Groenen, [2005] Borg, I. and Groenen, P. J.  _Modern Multidimensional Scaling: Theory and Applications_.  Springer Science & Business Media, New York, NY, 2005.
  * Duin and Pekalska, [2005] Duin, R. P. and Pekalska, E.  _The Dissimilarity Representation for Pattern Recognition: Foundations and Applications_ , volume 64.  World Scientific, Singapore, 2005.
  * Platt, [2005] Platt, J. C.  Fastmap, Metricmap, and Landmark MDS are all Nyström algorithms.  In _Proc. the 10th International Workshop on Artificial Intelligence and Statistics_ , pages 261–268. 2005.
  * Law and Jain, [2006] Law, M. and Jain, A.  Incremental nonlinear dimensionality reduction by manifold learning.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 28:377–391, 2006.
  * Takeuchi and Yamanishi, [2006] Takeuchi, J. and Yamanishi, K.  A unifying framework for detecting outliers and change points from time series.  _IEEE Transactions on Knowledge and Data Engineering_ , 18:482–492, 2006.
  * Adams and MacKay, [2007] Adams, R. P. and MacKay, D. J.  Bayesian online changepoint detection.  _Preprint arXiv:0710.3742_ , 2007.
  * Lee et al., [2007] Lee, J. A., Verleysen, M., et al.  _Nonlinear Dimensionality Reduction_ , volume 1.  Springer, New York, NY, 2007.
  * Van der Maaten and Hinton, [2008] Van der Maaten, L. and Hinton, G.  Visualizing data using t-SNE.  _Journal of Machine Learning Research_ , 9:2579–2605, 2008.
  * Schuon et al., [2008] Schuon, S., Durković, M., Diepold, K., Scheuerle, J., and Markward, S.  Truly incremental locally linear embedding.  In _Proc. CoTeSys 1st International Workshop on Cognition for Technical Systems_. 2008.
  * Wilson and Zhu, [2008] Wilson, R. C. and Zhu, P.  A study of graph spectra for comparing graphs and trees.  _Pattern Recognition_ , 41:2833–2841, 2008.
  * Chandola et al., [2009] Chandola, V., Banerjee, A., and Kumar, V.  Anomaly detection: A survey.  _ACM Computing Surveys_ , 41:1–58, 2009.
  * Van Der Maaten et al., [2009] Van Der Maaten, L., Postma, E., Van den Herik, J., et al.  Dimensionality reduction: a comparative review.  Technical report, Tilburg University, 2009.
  * Abdi and Williams, [2010] Abdi, H. and Williams, L. J.  Principal component analysis.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 2:433–459, 2010.
  * Durrant and Kabán, [2010] Durrant, R. J. and Kabán, A.  Compressed fisher linear discriminant analysis: Classification of randomly projected data.  In _Proc. the 16th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 1119–1128. 2010.
  * Barthélemy, [2011] Barthélemy, M.  Spatial networks.  _Physics Reports_ , 499:1–101, 2011.
  * Lopes et al., [2011] Lopes, M., Jacob, L., and Wainwright, M. J.  A more powerful two-sample test in high dimensions using random projection.  In _Proc. Advances in Neural Information Processing Systems_ , volume 24, pages 1206–1214. 2011.
  * Halmos, [2012] Halmos, P. R.  _Finite Dimensional Vector Spaces, Reprinted Version_.  Martino Publishing, Mansfield Centre, CT, 2012.
  * Buitinck et al., [2013] Buitinck, L., Louppe, G., Blondel, M., Pedregosa, F., Mueller, A., Grisel, O., Niculae, V., Prettenhofer, P., Gramfort, A., Grobler, J., Layton, R., VanderPlas, J., Joly, A., Holt, B., and Varoquaux, G.  API design for machine learning software: experiences from the scikit-learn project.  In _Proc. ECML PKDD Workshop: Languages for Data Mining and Machine Learning_ , pages 108–122. 2013.
  * Van Der Maaten, [2013] Van Der Maaten, L.  Barnes-hut-sne.  _Preprint arXiv:1301.3342_ , 2013.
  * Vanhems et al., [2013] Vanhems, P., Barrat, A., Cattuto, C., Pinton, J.-F., Khanafer, N., Régis, C., Kim, B.-a., Comte, B., and Voirin, N.  Estimating potential infection transmission routes in hospital wards using wearable proximity sensors.  _PLoS ONE_ , 8:e73970, 2013.
  * Wang et al., [2013] Wang, D., Qiao, H., Zhang, B., and Wang, M.  Online support vector machine based on convex hull vertices selection.  _IEEE Transactions on Neural Networks and Learning Systems_ , 24:593–609, 2013.
  * Yang et al., [2013] Yang, Z., Peltonen, J., and Kaski, S.  Scalable optimization of neighbor embedding for visualization.  In _Proc. the 30th International Conference on Machine Learning_ , pages 127–135. 2013.
  * Hoi et al., [2014] Hoi, S. C., Wang, J., and Zhao, P.  Libol: A library for online learning algorithms.  _Journal of Machine Learning Research_ , 15:495–499, 2014.
  * Akoglu et al., [2015] Akoglu, L., Tong, H., and Koutra, D.  Graph based anomaly detection and description: a survey.  _Data Mining and Knowledge Discovery_ , 29:626–688, 2015.
  * Boutsidis et al., [2015] Boutsidis, C., Garber, D., Karnin, Z., and Liberty, E.  Online principal components analysis.  In _Proc. the 2015 Annual ACM-SIAM Symposium on Discrete Algorithms_ , pages 887–901. 2015.
  * Cunningham and Ghahramani, [2015] Cunningham, J. P. and Ghahramani, Z.  Linear dimensionality reduction: survey, insights, and generalizations.  _Journal of Machine Learning Research_ , 16:2859–2900, 2015.
  * Gisbrecht and Hammer, [2015] Gisbrecht, A. and Hammer, B.  Data visualization by nonlinear dimensionality reduction.  _WIREs Data Mining and Knowledge Discovery_ , 5:51–73, 2015.
  * Gisbrecht and Schleif, [2015] Gisbrecht, A. and Schleif, F.-M.  Metric and non-metric proximity transformations at linear costs.  _Neurocomputing_ , 167:643–657, 2015.
  * Gisbrecht et al., [2015] Gisbrecht, A., Schulz, A., and Hammer, B.  Parametric nonlinear dimensionality reduction using kernel t-SNE.  _Neurocomputing_ , 147:71–82, 2015.
  * Juvonen et al., [2015] Juvonen, A., Sipola, T., and Hämäläinen, T.  Online anomaly detection using dimensionality reduction techniques for HTTP log analysis.  _Computer Networks_ , 91:46–56, 2015.
  * Schleif and Tino, [2015] Schleif, F.-M. and Tino, P.  Indefinite proximity learning: A review.  _Neural Computation_ , 27:2039–2096, 2015.
  * Gepperth and Hammer, [2016] Gepperth, A. and Hammer, B.  Incremental learning algorithms and applications.  In _Proc. European Symposium on Artificial Neural Networks_ , pages 357–368. 2016.
  * Liberty et al., [2016] Liberty, E., Sriharsha, R., and Sviridenko, M.  An algorithm for online k-means clustering.  In _Proc. the Meeting on Algorithm Engineering and Experiments_, pages 81–89. 2016.
  * Malik et al., [2016] Malik, Z. K., Hussain, A., and Wu, J.  An online generalized eigenvalue version of laplacian eigenmaps for visual big data.  _Neurocomputing_ , 173:127–136, 2016.
  * Zhou et al., [2016] Zhou, X., Zhang, X., and Wang, B.  Online support vector machine: A survey.  In _Proc. Harmony Search Algorithm: the 2nd International Conference on Harmony Search Algorithm_ , pages 269–278. 2016.
  * Einsiedler et al., [2017] Einsiedler, M., Ward, T., et al.  _Functional Analysis, Spectral Theory, and Applications_ , volume 276.  Springer, Cham, Switzerland, 2017.
  * Rahman, [2017] Rahman, M. S.  _Basic Graph Theory_.  Springer, Cham, Switzerland, 2017.
  * Cardot and Degras, [2018] Cardot, H. and Degras, D.  Online principal component analysis in high dimension: Which algorithm to choose?  _International Statistical Review_ , 86:29–50, 2018.
  * Donnat and Holmes, [2018] Donnat, C. and Holmes, S.  Tracking network dynamics: A survey using graph distances.  _Annals of Applied Statistics_ , 12:971–1012, 2018.
  * P. Peixoto and Gauvin, [2018] P. Peixoto, T. and Gauvin, L.  Change points, memory and epidemic spreading in temporal networks.  _Scientific Reports_ , 8:15511, 2018.
  * Hjuler et al., [2019] Hjuler, N., Italiano, G. F., Parotsidis, N., and Saulpic, D.  Dominating sets and connected dominating sets in dynamic graphs.  In _Proc. the 36th International Symposium on Theoretical Aspects of Computer Science_ , volume 126, pages 35:1–35:17. 2019\.
  * Long and Ferguson, [2019] Long, A. W. and Ferguson, A. L.  Landmark diffusion maps (L-dMaps): Accelerated manifold learning out-of-sample extension.  _Applied and Computational Harmonic Analysis_ , 47:190–211, 2019\.
  * Masuda and Holme, [2019] Masuda, N. and Holme, P.  Detecting sequences of system states in temporal networks.  _Scientific Reports_ , 9:795, 2019.
  * Velliangiri et al., [2019] Velliangiri, S., Alagumuthukrishnan, S., and Thankumar joseph, S. I.  A review of dimensionality reduction techniques for efficient computation.  _Procedia Computer Science_ , 165:104–111, 2019.
  * Ayesha et al., [2020] Ayesha, S., Hanif, M. K., and Talib, R.  Overview and comparative study of dimensionality reduction techniques for high dimensional data.  _Information Fusion_ , 59:44–58, 2020.
  * Masuda and Lambiotte, [2020] Masuda, N. and Lambiotte, R.  _A Guide to Temporal Networks_.  World Scientific, Singapore, 2nd edition, 2020.
  * Porter, [2020] Porter, M. A.  _Emerging Frontiers in Nonlinear Science_ , pages 131–159.  Springer, Cham, Switzerland, 2020.
  * Ahmad et al., [2021] Ahmad, W., Porter, M. A., and Beguerisse-Díaz, M.  Tie-decay networks in continuous time and eigenvector-based centralities.  _IEEE Transactions on Network Science and Engineering_ , 8:1759–1771, 2021.
  * Böckenhauer et al., [2021] Böckenhauer, H.-J., Hromkovič, J., Krug, S., and Unger, W.  On the advice complexity of the online dominating set problem.  _Theoretical Computer Science_ , 862:81–96, 2021.
  * Hoi et al., [2021] Hoi, S. C., Sahoo, D., Lu, J., and Zhao, P.  Online learning: A comprehensive survey.  _Neurocomputing_ , 459:249–289, 2021.
  * Migenda et al., [2021] Migenda, N., Möller, R., and Schenck, W.  Adaptive dimensionality reduction for neural network-based online principal component analysis.  _PLoS ONE_ , 16:e0248896, 2021.
  * Sugishita et al., [2021] Sugishita, K., Porter, M. A., Beguerisse-Díaz, M., and Masuda, N.  Opinion dynamics on tie-decay networks.  _Physical Review Research_ , 3:023249, 2021.
  * Zhang et al., [2021] Zhang, H., Wang, P., Gao, X., Qi, Y., and Gao, H.  Out-of-sample data visualization using bi-kernel t-SNE.  _Information Visualization_ , 20:20–34, 2021.
  * Zuo and Porter, [2021] Zuo, X. and Porter, M. A.  Models of continuous-time networks with tie decay, diffusion, and convection.  _Physical Review E_ , 103:022304, 2021.
  * Ghojogh et al., [2023] Ghojogh, B., Crowley, M., Karray, F., and Ghodsi, A.  _Elements of Dimensionality Reduction and Manifold Learning_.  Springer Nature, Cham, Switzerland, 2023.
  * Thongprayoon et al., [2023] Thongprayoon, C., Livi, L., and Masuda, N.  Embedding and trajectories of temporal networks.  _IEEE Access_ , 11:41426–41443, 2023.