Online Sparse Subspace Clustering

TL;DR

This paper introduces an online algorithm for sparse subspace clustering that efficiently handles streaming data by updating similarities and clustering results without revisiting the entire dataset.

Contribution

It proposes a novel online approach for sparse subspace clustering with convergence guarantees and dynamic regret analysis for non-strongly convex costs.

Findings

Online algorithm converges to a neighborhood of batch solution

Dynamic regret bounds are established for non-strongly convex costs

Method enables real-time clustering of streaming data

Abstract

This paper focuses on the sparse subspace clustering problem, and develops an online algorithmic solution to cluster data points on-the-fly, without revisiting the whole dataset. The strategy involves an online solution of a sparse representation (SR) problem to build a (sparse) dictionary of similarities where points in the same subspace are considered "similar," followed by a spectral clustering based on the obtained similarity matrix. When the SR cost is strongly convex, the online solution converges to within a neighborhood of the optimal time-varying batch solution. A dynamic regret analysis is performed when the SR cost is not strongly convex.

Equations

$$C_t^{\ast }\in \underset{C\in {\mathbb{R}}^{T\times T}}{\mathrm{argmin}}$$

$$\text{diag}(C)=0$$

$$C_{n+1}={\text{prox}}_{\gamma \Vert \cdot {\Vert }_1,\text{diag}(\cdot )=0}\circ (I-\gamma \nabla f_t)(C_n)$$

$$C_{t+1}=\underset{P_t}{\underbrace{{\text{prox}}_{{\gamma }_t\Vert \cdot {\Vert }_1,\text{diag}(\cdot )=0}}}\circ \underset{G_t}{\underbrace{(I-{\gamma }_t\nabla f_t)}}(C_t)\text{vspace}-.9em$$

$$\forall t\ge 1,\Vert C_t-C_t^{\ast }\Vert \le {\tilde{L}}_{t-1}\left(\Vert C_0-C_0^{\ast }\Vert +\sum _{\tau =0}^{t-1}\frac{{\sigma }_{\tau }}{{\tilde{L}}_{\tau }}\right)$$

$$\Vert G$$

$$=\Vert C-{\gamma }_t\nabla f_t(C)-C^{\prime }+{\gamma }_t\nabla f_t(C^{\prime }){\Vert }^2$$

$$\le (1-2{\gamma }_t{m}_t+{\gamma }_t^2{M}_t^2)\Vert C-C^{\prime }{\Vert }^2$$

$$\le L_t^2\Vert C-C^{\prime }{\Vert }^2.$$

$$\Vert C_{t+1}-C_{t+1}^{\ast }\Vert$$

$$\le L_t\Vert C_t-C_t^{\ast }\Vert +{\sigma }_t$$

$$\Vert C_t-C_t^{\ast }\Vert \le {\left({\hat{L}}_{t-1}\right)}^t\Vert C_0-C_0^{\ast }\Vert +\frac{{\hat{\sigma }}_t}{1-{\hat{L}}_{t-1}}.$$

$$\displaystyle\frac{1}{t}\sum_{\tau=0}^{t-1}\left(F_{\tau}(C_{\tau})-F_{\tau}(C_{\tau}^{*})\right)\leq\frac{M}{2t}\|C_{0}-\hat{C}_{t}^{*}\|^{2}+\frac{1}{t}\Biggl{(}F_{0}(C_{0})$$

$$\displaystyle\hskip 28.45274pt-F_{t-1}(C_{t})+\sum_{{\tau}=0}^{t-2}\delta_{\tau}+\sum_{{\tau}=0}^{t-1}\rho_{t}({\tau})(T+\frac{M}{2}\rho_{t}({\tau}))\Biggr{)}$$

$$\frac{1}{t}\sum _{\tau =0}^{t-1}\left(F_{\tau }(C_{\tau +1})-F_{\tau }({\hat{C}}_t^{\ast })\right)\le \frac{M}{2t}\Vert C_0-{\hat{C}}_t^{\ast }{\Vert }^2$$

$$\frac{1}{t}\sum _{\tau =0}^{t-1}\left(F_{\tau }({\hat{C}}_t^{\ast })-F_{\tau }(C_{\tau }^{\ast })\right)\le \frac{1}{t}\sum _{\tau =0}^{t-1}{\rho }_t(\tau )(T+\frac{M}{2}{\rho }_t(\tau ))$$

$$\frac{1}{t}\sum _{\tau =0}^{t-1}\left(F_{\tau }(C_{\tau })-F_{\tau }(C_{\tau +1})\right)=\frac{1}{t}(F_0(C_0)-F_{t-1}(C_t)$$

$$+\sum _{\tau =0}^{t-2}\left(f_{\tau +1}(C_{\tau +1})-f_{\tau }(C_{\tau +1})\right)$$

$$\le \frac{1}{t}\left(F_0(C_0)-F_{t-1}(C_t)+\sum _{\tau =0}^{t-2}{\delta }_{\tau }\right)$$

$$\displaystyle\frac{1}{t}\sum_{{\tau}=0}^{t-1}\left(F_{\tau}(C_{\tau})-F_{\tau}(C_{\tau}^{*})\right)\leq\frac{1}{t}\biggl{(}\frac{M}{2}\|C_{0}-\hat{C}_{t}^{*}\|^{2}$$

$$\displaystyle+F_{0}(C_{0})-F_{t-1}(C_{t})\biggr{)}+\frac{M\hat{\rho}_{t}^{2}}{2}+T\hat{\rho}_{t}+\hat{\delta}_{t}.$$

Full text

**Index Terms— ** Subspace clustering; sparse representation; time-varying optimization; online algorithms.

1 Introduction and Problem Setup

Modern data processing tasks aim to extract information from datasets or signals on graphs – examples include identification of trends or patterns, learning of dynamics and data structures, or methods for comprehensive awareness of the underlying networks or systems generating the data [1, 2]. In this domain, the present paper focuses on data processing methods for streams of (possibly high-dimensional) data, with particular emphasis on a setting where underlying computational constraints require one to process data on-the-fly, and with limited access to stored data.

One prominent task is clustering, which is utilized to cluster data points based on well-defined metrics modeling similarities (or distances) among data points, or capturing underlying data structures. For example, spectral clustering groups data points based on minimizing cuts of a similarity graph [3]. In particular, subspace clustering builds a similarity graph where points in the same subspace are considered “similar”. It does this by finding either a low-rank representation (LRR) [4, 5] or sparse representation (SSC) [6, 7] of the data such that the data points are represented as linear combinations of other data points in the same subspace. In this paper, we will consider an online SSC methodology, where underlying computational considerations prevent one from solving pertinent optimization problems associated with a given set of data [6] before a new datum arrives.

To outline the problem concretely, consider $N$-dimensional data points $\{x_{t},t\in\mathbb{N}\}$ sequentially arriving at times $t\cdot h$, with $h>0$ a given interval111Notation: hereafter, $(\cdot)^{T}$ denotes transposition. For a given vector $x\in\mathbb{R}^{N}$ or matrix $X\in\mathbb{R}^{N\times M}$, $\|x\|$ and $\|X\|$ refer to a generic norm, and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|X\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ the spectral norm. If $X_{ij}$ is the $(i,j)$ entry of $X$, then $\|X\|_{F}^{2}=\sum_{i,j}X_{ij}^{2}$ and $\|X\|_{1}=\sum_{i,j}|X_{ij}|$. The composition of two operators is denoted by $\circ$. Write $f(n)=\mathcal{O}(g(n))$ to denote that for sufficiently large $n$, $\exists c>0$ such that $|f(n)|\leq c|g(n)|$. . Assume that the data points lie in (or in the neighbourhood of) $S$ subspaces $\{\mathcal{S}_{i}\}_{i=1}^{S}\subset\mathbb{R}^{N}$, with $\text{dim}(\mathcal{S}_{i})=d_{i}$ for each $i$. This paper studies the problem of repeatedly applying SSC to all observed data up to a given time, $\{x_{j},j\leq t\}$ with $\bar{X}_{t}=[x_{1},\ldots,x_{t}]$. SSC is a two-step approach [6]: first, step [S1], based on the self-expressiveness property, a sparse representation (SR) problem is solved to identify the (sparse) coefficients $\{\bar{c}_{j}\}_{j=1}^{t}$ so that $x_{j}=\bar{X}_{t}\bar{c}_{j}$ for all $j=1,\ldots,t$; that is, data points are represented as linear combinations of data points in the same subspace (and we force the $j^{\text{th}}$ component of $\bar{c}_{j}$ to be [math] to exclude the trivial solution). Second, step [S2], apply spectral clustering based on a similarity matrix $W:=|\bar{C}_{t}|+|\bar{C}_{t}^{T}|$, where $\bar{C}_{t}:=[\bar{c}_{1},\ldots\bar{c}_{t}]$ and $|\cdot|$ is taken entry-wise. However, in the streaming setting, this setting has two drawbacks:

1. (d2)

The dimensions of $\bar{X}_{t}$ and $\bar{C}_{t}$ grow with time, thus increasing the complexity of the associated SR and spectral clustering tasks; and, 2. (d2)

Due to underlying computational complexity considerations, steps [S1] and [S2] might not be executed to completion within a time interval $h$ (i.e., before a new datum arrives).

Given (d1)-(d2), we address the problem of developing online algorithmic solutions to carry out steps [S1]-[S2] at each time $t$, based on a given computational budget. The first step towards this goal involves the processing of data points using a “sliding window” $X_{t}:=[x_{t-T+1},\ldots,x_{t}]\in\mathbb{R}^{N\times T}$ of length $T$ (with $T$ determined by the computational budget, as explained later in the paper). Step [S1] is ideally carried out at time $t$ by solving the following SR problem:

[TABLE]

with $\lambda_{e}>0$ a given tuning parameter. If $\lambda_{e}$ is too small, in particular if $\lambda_{e}\leq\|\textrm{vec}(X_{t}^{T}X_{t})\|_{\infty}^{-1}$, then the minimizer $C_{t}^{*}$ may have all-zero columns, which is not informative, hence we always choose $\lambda_{e}$ sufficiently large.

Solving (SRt) to convergence within a time interval $h$ might not be possible for a given computational budget, especially for streams of high-dimensional vectors over a large window $T$. Section 2 will address the design of online algorithms, for the case where only one algorithmic step can be performed before a new datum $x_{t}$ arrives (the case of multiple steps follows easily).

With the minimizer (or approximate solution) $C_{t}^{*}$, one can compute the matrix $W_{t}=|C_{t}^{*}|+|C_{t}^{*}|^{T}$. Interpreting $W$ as the “similarity” matrix of a graph, one can compute the graph Laplacian $L_{t}=D_{t}-W_{t}$ where the degree matrix $D_{t}$ is a diagonal matrix attained by summing the rows of $W_{t}$. The graph Laplacian is then normalized in one of two possible ways: as the symmetric graph Laplacian, $L_{sym,t}=D_{t}^{-1/2}L_{t}D_{t}^{-1/2}$, or as the random walk graph Laplacian, $L_{rw,t}=D_{t}^{-1}L_{t}$. We then compute the $S$ trailing eigenvectors of the normalized graph Laplacian and, viewing these eigenvectors as columns, we cluster their rows in $\mathbb{R}^{S}$ into $S$ cluster using the k-means algorithm; see, e.g., [3] for details on spectral clustering. This clustering is then applied to the original data points. Section 3 will elaborate further on this step.

2 Online Sparse Representation

The proximal gradient descent algorithm and its accelerated version [8] have rigorous convergence guarantees and can be applied to equations of the form (SRt), which we now detail. Let $f_{t}(C)=\frac{\lambda_{e}}{2}\|X_{t}-X_{t}C\|^{2}_{F}$ for brevity, and notice that $\nabla f_{t}(C)=\lambda_{e}X_{t}^{T}(X_{t}C-X_{t})$ so $\nabla f_{t}$ is Lipschitz continuous with constant $M_{t}=\lambda_{e}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|X_{t}^{T}X_{t}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$. Let $n$ be the iteration index of the algorithm, and let $\gamma<\frac{2}{M_{t}}$. Then the (batch) proximal gradient descent algorithm, used to solve (SRt) at time $t$, involves the following iterations for $n=1,2,\dots$ until convergence:

[TABLE]

where $\text{prox}_{g,\mathcal{X}}(z)=\text{argmin}_{x\in\mathcal{X}}g(x)+\frac{1}{2}\|x-z\|^{2}$ is the proximal operator defined over a closed convex set $\mathcal{X}$ and for a function $g$. Convergence of (1) to a minimizer $C_{t}^{*}$ is shown in, e.g., [9]. Furthermore, $F_{t}(C_{n})-F_{t}(C_{t}^{*})\rightarrow 0$ by the continuity of $F_{t}$; see Theorem 10.21 in [8] for the rate.

Consider now the case where only one iteration (1) can be performed per time interval $h$ (see Remark 3 for the extension to multiple steps). Then, an online implementation of the proximal gradient descent algorithm involves the sequential execution of the following step at each time $t$:

[TABLE]

where the coefficient $\gamma_{t}$ is selected so that $\gamma_{t}<\frac{2}{M_{t}}$. The difference between (1) and (2) is that $f_{t}$ changes per iteration in the latter. The goal is to demonstrate that the online algorithm (2) can track the sequence of optimizers $\{C_{t}^{*}\}$. In the following, the performance of the online algorithm (2) is investigated in two cases:

i) The cost of (SRt) is strongly convex for each $t$. In this case, we will derive bounds for $\|C_{t}-C_{t}^{*}\|$, where $C_{t}^{*}$ is the unique minimizer of (SRt) for each $t$.

ii) The cost of (SRt) is not strongly convex. In this case, we will derive dynamic regret bounds.

Before proceeding, to capture the variability of the clustering solutions, define $\sigma_{t}=\|C_{t+1}^{*}-C_{t}^{*}\|_{F}$ [10]. As expected, it will be shown that high variability leads to poor tracking performance. The following assumptions are then introduced.

Assumption 1**.**

The matrix $X_{t}^{T}X_{t}$ is positive definite for each time $t$ (e.g., $N>T$ and $X_{t}$ is full rank). Let $m_{t}>0$ be the smallest eigenvalue of $X_{t}^{T}X_{t}$.

By inspection, the $m_{t}$ of Assumption 1 is the strong convexity constant of $f_{t}$. Based on the assumptions below, the first result is stated next, where $\|\cdot\|$ is taken to be the Frobenius norm.

Theorem 1**.**

Under Assumption 1,

[TABLE]

where $L_{t}=\max\{|1-\gamma_{t}m_{t}|,|1-\gamma_{t}M_{t}|\}$, $\tilde{L}_{t}=\prod_{\tau=0}^{t}L_{\tau}$.

Proof.

Define $G_{t}$ and $P_{t}$ as in (2), so $C_{t+1}=(P_{t}\circ G_{t})C_{t}$. Observe, for any $C,C^{\prime}\in\mathbb{R}^{T\times T}$,

[TABLE]

Also, by optimality, $C_{t}^{*}$ is a fixed point, so $C_{t}^{*}=P_{t}G_{t}C_{t}^{*}$. Therefore, one has that $\|C_{t+1}-C_{t}^{*}\|=\|P_{t}G_{t}C_{t}-P_{t}G_{t}C_{t}^{*}\|\leq\|G_{t}C_{t}-G_{t}C_{t}^{*}\|\leq L_{t}\|C_{t}-C_{t}^{*}\|$, where the second inequality comes from the nonexpansiveness of the prox operator [9]. Finally, we get

[TABLE]

and we apply this inequality recursively. ∎

Corollary 1**.**

Define $\hat{L}_{t}=\underset{\tau=0,...,t}{\max}L_{\tau}$ and $\hat{\sigma}_{t}=\underset{\tau=0,...,t}{\max}\sigma_{\tau}$. If Assumption 1 holds, then, for each $t$:

[TABLE]

If $m_{\tau}\geq m$ and $M_{\tau}\leq M$ and $\gamma_{\tau}$ is chosen $\gamma_{\tau}=2/(m_{\tau}+M_{\tau})$ for all $\tau=0,\ldots,t$, then $\hat{L}_{t}\leq(M-m)/(M+m)<1$.

Proof.

Follows from Thm. 1 and the geometric series. ∎

These results closely follow the analysis of [10, 11, 12], applying it to SSC.

Remark 1**.**

Under Assumption 1, it can be shown that $F_{t}(C_{t})-F_{t}(C_{t}^{*})\leq\frac{M_{t}}{2}\|C_{t}-C_{t}^{*}\|^{2}$ [8, Thm. 10.29].

In the following, we consider the case where the cost of (SRt) is not strongly convex. It is clear that contractive arguments cannot be utilized in this case since (2) is no longer a strongly monotone operator. Again, we use $\|\cdot\|$ for the Frobenius norm, and $\|g\|_{\infty}=\sup_{x}\,|g(x)|$.

Theorem 2**.**

Let $\hat{C}_{t}^{*}=\frac{1}{t}\sum_{\tau=0}^{t-1}C_{\tau}^{*}$, $\rho_{t}(\tau)=\|\hat{C}_{t}^{*}-C_{\tau}^{*}\|$, $\delta_{t}=\|f_{t+1}-f_{t}\|_{\infty}$, and $M=\max_{t}M_{t}$. Set $\gamma_{t}=\frac{1}{M}$. Then

[TABLE]

Proof.

From the descent lemma [8, Thm. 10.16], we have

[TABLE]

Using the Lipschitz continuity of $\nabla f_{t}$, $\partial\|\cdot\|_{1}\in[-1,1]^{T^{2}}$, and the Cauchy-Schwartz inequality, we get

[TABLE]

And, rearranging, we find

[TABLE]

Adding these together gives us the result. ∎

Corollary 2**.**

*Define $\hat{\rho}_{t}=\underset{\tau=0,...,t-1}{\max}\rho_{t}(\tau)$ and

$\hat{\delta}_{t}=\underset{\tau=0,...,t-2}{\max}\delta_{\tau}$. Then*

[TABLE]

The term $\hat{C}_{t}^{*}$ serves as a center point of all the $C^{*}_{\tau}$ and is the most meaningful “best overall” point. The corollary says that if $\hat{\delta}_{t}$ and $\hat{\rho}_{t}$ are well-behaved (i.e., the function changes slowly) then, on average, $F_{t}(C_{t})$ tracks within a constant term of $F_{t}(C_{t}^{*})$.

Remark 2**.**

If $(C_{t}^{*})$ is bounded, then $\hat{\rho}_{t}$ converges and $\hat{C}_{t}^{*}$ has a convergent subsequence. In that case, we can replace them with their limits in the bound in Corollary 2. The bound is only meaningful, though, when the $\delta_{t}$’s are finite. One way to make them finite is to impose boundedness with respect to the infinity norm as another constraint in (SRt). This particular constraint can be incorporated into our closed-form proximal projection [9, Prop. 24.47].

Remark 3**.**

If $f_{t}$ is not strongly convex (which it will not be if $T>N$), then we can add a Tikhonov term to make it strongly convex. For example, we could add $\frac{\lambda_{r}}{2}\|C\|^{2}$. While this will provide the stronger result of Theorem 1, it incurs an error in the minimizer. That is, $\|C_{r,t}-C_{t}^{*}\|\leq\|C_{r,t}-C_{r,t}^{*}\|+\|C_{r,t}^{*}-C_{t}^{*}\|$ where $C_{r,t}$ and $C_{r,t}^{*}$ are the regularized tracking sequence and regularized minimizer sequence respectively.

Remark 4**.**

If we take more than one iteration per time step, then we can modify $(f_{t})$ accordingly in order to use the results in this paper. For example, with 2 iterations per time step, define $\tilde{f}_{t}=f_{\left\lfloor t/2\right\rfloor}$. Alternatively, we can modify the Theorems. For example, if we take $n_{t}$ iterations at time $t$, then we just have to redefine $\tilde{L}_{t}=\prod_{\tau=0}^{t}L_{\tau}^{n_{\tau}}$ in Theorem 1.

Note that we did not consider accelerating our algorithm. In the non-strongly convex case, accelerated methods rely on global structure not just local descent. Because of this, it is not obvious that adapting an accelerated method to the online setting would lead to better tracking. However, this is a further research direction that we are currently exploring.

3 Spectral Clustering

The main factor in determining how many iterations to take in step [S1] is the ratio between the costs of steps [S1] and [S2]. There are papers that explore online spectral clustering [14], but the results require relatively small changes in the graph. There are no such guarantees here. For example, when the oldest data point is thrown out and a new one added, all connections to the previous data point are thrown out as well, new connections have to be made for all of the points that were previously connected to the old data point, and connections have to be made for the new data point. While the end result of the spectral clustering may barely change, the change in the graph is catastrophic. Thus, we do a full batch spectral clustering operation, and leave a more general framework for online spectral clustering as a future research direction.

Note that the proximal operator in (2) simplifies to $\text{proj}_{\text{diag}(\cdot)=0}\circ\text{prox}_{\gamma_{t}\|\cdot\|_{1}}$. The closed form expression for the latter proximal operator is soft-thresholding: $\text{prox}_{\gamma_{t}\|\cdot\|_{1}}C=\text{sign}(C)(|C|-\gamma)_{+}$ component-wise where $(a)_{+}=\text{max}(a,0)$. Thus, the cost of each iteration of step [S1] is dominated by the gradient descent sub-step, and so the total cost of each step is $\mathcal{O}(NT^{2})$ operations. If, instead of applying the proximal gradient operator just once to $C_{t}$, we apply it $n_{t}$ times, then the step at time $t$ will cost $\mathcal{O}(n_{t}NT^{2})$ operations.

The cost of step [S2] depends on what method we use to compute the specified eigenvectors. An upper bound on the cost is $\mathcal{O}(T^{3})$ which can be achieved by classical dense algorithms and gives all $T$ eigenvalues $\lambda$. To find $S\ll T$ eigenvalues approximately, the power method or Lanczos iterative methods can be used [15]. First consider that $L_{rw}v=\lambda v$ iff $v-D^{-1}Wv=\lambda v$ iff $D^{-1}Wv=(1-\lambda)v$. Also, by the Gershgorin disc theorem, the eigenvalues of $L_{rw}$ are in [0,2] and the eigenvalues of $D^{-1}W$ are in [-1,1]. Thus, for $D^{-1}W$, we want the eigenvectors corresponding to the eigenvalues near 1. The leading eigenpairs of $D^{-1}W$ should correspond to the trailing eigenpairs of $L_{rw}$ as long as the trailing eigenvalues of $L_{rw}$ are closer to 0 than the leading eigenvalues are to 1. If this is not the case, then we have to use the power method to compute more than S eigenpairs, so that we can take the eigenvectors corresponding to the S largest positive eigenvalues. The power method costs $\mathcal{O}(nnz\cdot niter\cdot S)$ in the ideal case where we do not have to compute extra eigenpairs. Here, $nnz$ denotes the number of nonzero elements and $niter$ denotes the number of iterations to reach convergence.

There are efficient methods for computing the trailing eigenpairs of $L_{rw}$ directly. In particular, [16] found that the Jacobi-Davidson method was superior to the Lanczos method, in terms of computation time, for spectral clustering. Both methods cost $\mathcal{O}(nnz\cdot niter)$, but $niter\approx S$ can vary.

The number of nonzeros, in the ideal case, can be estimated by the subspace dimensions, $d_{i}$. For a given subspace, the minimum number of points needed to represent another point, as long as it isn’t cohyperplanar with a strict subset of the points, is exactly the dimension of the subspace. Thus, we can say that $nnz(C^{*})\geq\mathcal{O}(\sum_{i}d_{i}T_{i}^{t})$ where $T_{i}^{t}$ is the number of data points in $\mathcal{S}_{i}$ at time $t$. The same can be said of $W$, $L$, and $L_{rw}$. In the case where each subspace has the same dimension, $d$, and the same number of points in it, $T/S$, then $nnz(C^{*})\geq\mathcal{O}(dT)$ (and recall $d\leq N$). This means in the best case, [S2] costs $\mathcal{O}(dTS)$, while a single step of [S1] costs $\mathcal{O}(NT^{2})$ so this suggests choosing $T\approx dS/N$ to balance the costs of the two steps; if $T$ is smaller than this, multiple steps to solve [S1] can be taken.

4 Numerical Results

We performed tests on both synthetic data and the Yale Face Database [17]. The synthetic data was composed of $S=10$ 5-dimensional subspaces in $R^{50}$ ($d=5$), each with 50 points. Noise was added to make the simulation more realistic. The sliding window had a capacity of $T=400$ data points and we took 100 time steps. The Yale Face Database has $S=38$ subspaces, and we let the sliding window have a capacity of $T=500$ points and took 200 time steps. The points in the Yale Face Database are in $\mathbb{R}^{2016}$. For both datasets, we took 50 iterations of the optimization algorithm per time step (cf. Remark 4).

For the synthetic data, $T>N$, so the cost function is not strongly convex. On the other hand, for the Yale Face Database, $T<N$, so the cost function is strongly convex. The numerical results show that for both datasets, though, the objective value trajectory converges to a region above the minimum trajectory. This can be seen in Figure 1.

The objective error seems to be driving whether or not the clustering error converges. Figure 2 shows the clustering error.

For both datasets, the clustering error of the tracking algorithm decreases. However, decreasing the number of iterations per time step too much causes the clustering error to no longer decrease. This suggests that there is some maximum value of the objective error for the clustering error to decrease.

Finally, for the synthetic dataset, step [S2] took 50 times as long as step [S1]. For the Yale Face dataset, step [S2] took 3 times as long as step [S1]. While we took a sufficient number of iterations in order for the clustering error of our algorithm to be small, for an actual system, its dynamics would dictate the number of iterations. The spectral clustering time would be subtracted from the length of a time step, and this value would be divided by the time for step [S1].