\section{Introduction}\label{sec:intro}
In {the} past a few years, subspace clustering~\cite{vidal2010tutorial,soltanolkotabi2012geometric} has been extensively studied and has established solid applications, for example, in computer vision~\cite{elhamifar2009sparse} and network topology inference~\cite{eriksson2011high}.  Among many subspace clustering algorithms which aim to obtain a structured representation  to fit the underlying data, two prominent examples are  Sparse Subspace Clustering~(SSC)~\cite{elhamifar2009sparse,soltanolkotabi2014robust} and Low-Rank Representation~(LRR)~\cite{liu2013robust}. Both of them utilize the idea of self-expressiveness, i.e., expressing each sample as a linear combination of the remaining. What is of difference is that SSC pursues a sparse solution while LRR prefers a low-rank structure.

In this paper, we are interested in the LRR method, which is shown to achieve state-of-the-art performance on a broad range of real-world problems~\cite{liu2013robust}. Recently, \cite{liu2014recovery} demonstrated that, when equipped with a proper dictionary, LRR can even handle the coherent data~--~a challenging issue in the literature~\cite{candes2009exact,candes2011robust} but commonly emerges in realistic datasets such as the Netflix.

Formally, the LRR problem we investigate here is formulated as follows~\cite{liu2013robust}:
\begin{align}
	\label{eq:lrr}
	\min_{\boldsymbol{X}, \boldsymbol{E}}\ \fractwo{\lambda_1} \fronorm{\boldsymbol{Z}-\boldsymbol{Y}\boldsymbol{X}-\boldsymbol{E}}^2 + \nuclearnorm{\boldsymbol{X}} + \lambda_2 \onenorm{\boldsymbol{E}}.
\end{align}
Here, $\boldsymbol{Z} = (\boldsymbol{z}_1, \boldsymbol{z}_2, \cdots, \boldsymbol{z}_n) \in \mathbb{R}^{p\times n}$ is the observation matrix with $n$ samples lying in a $p$-dimensional subspace. The matrix $\boldsymbol{Y} \in \mathbb{R}^{p\times n}$ is a given dictionary, $\boldsymbol{E}$ is some possible sparse corruption and $\lambda_1$ and $\lambda_2$ are two tunable parameters. Typically, $\boldsymbol{Y}$ is chosen as the dataset $\boldsymbol{Z}$ itself. The program seeks a low-rank representation $\boldsymbol{X} \in \mathbb{R}^{n\times n}$ among all samples, each of which can be approximated by a linear combination of the atoms in the dictionary $\boldsymbol{Y}$.

While LLR is mathematically elegant, three issues are immediately incurred in the face of big data:
\begin{issue}[Memory cost of $\boldsymbol{X}$]\label{is:memX}
	In the LRR formulation~\eqref{eq:lrr}, there is typically no sparsity assumption on $\boldsymbol{X}$. Hence, the memory footprint of $\boldsymbol{X}$ is proportional to $n^2$ which precludes most of the recently developed nuclear norm solvers~\cite{lin2010augmented,jaggi2010simple,avron2012efficient,hsieh2014nuclear}.
\end{issue}
\begin{issue}[Computational cost of $\nuclearnorm{\boldsymbol{X}}$]\label{is:compX}
	Due to the size of the nuclear norm regularized matrix $\boldsymbol{X}$ is $n \times n$, optimizing such problems can be computationally expensive even when $n$ is not too large~\cite{recht2010guaranteed}.
\end{issue}
\begin{issue}[Memory cost of $\boldsymbol{Y}$]\label{is:memY}
	Since  the dictionary size is $p\times n$, it is prohibitive to store the entire dictionary $\boldsymbol{Y}$ during optimization when  manipulating a huge volume of data.
\end{issue}
To remedy these issues, especially the memory bottleneck, one potential way is solving the problem in online manner. That is, we sequentially reveal the samples $\boldsymbol{z}_1, \boldsymbol{z}_2, \cdots, \boldsymbol{z}_n$ and update the components in $\boldsymbol{X}$ and $\boldsymbol{E}$. Nevertheless, such strategy appears difficult to execute due the the residual term in~\eqref{eq:lrr}. To be more precise, we note that each column of $\boldsymbol{X}$ is the coefficients of a sample with respect to the {\em entire} dictionary $\boldsymbol{Y}$, e.g., $\boldsymbol{z}_1 \approx \boldsymbol{Y} \boldsymbol{x}_1 + \boldsymbol{e}_1$. This indicates that without further technique, we have to load the entire dictionary $Y$ so as to update the columns of $\boldsymbol{X}$. Hence, for our purpose, we need to tackle a {more serious} challenge:
\begin{issue}[Partial realization of $\boldsymbol{Y}$]\label{is:partialY}
	We are required to guarantee the optimality of the solution but can only access part of the atoms of $Y$ in each iteration.
\end{issue}

\subsection{Related Works}
There are a vast body of works attempting to mitigate the memory and computational bottleneck of the nuclear norm regularizer. However, to the best of our knowledge, none of them can handle Issue~\ref{is:memY} and Issue~\ref{is:partialY} in the LRR problem.

One of the most popular ways to alleviate the huge memory cost is online implementation. \cite{feng2013online} devised an online algorithm for the Robust Principal Component Analysis~(RPCA) problem, which makes the memory cost independent of the sample size. Yet, compared to RPCA where the size of the nuclear norm regularized matrix is $p \times n$, that of LRR is $n\times n$~--~a worse and more challenging case. Moreover, their algorithm cannot address the partial dictionary issue that emerges in our case. It is also worth mentioning that~\cite{qiu2014recursive} established another online variant of RPCA. But since we are dealing with a different problem setting, i.e., the multiple subspaces regime, it is not clear how to extend their method to LRR.

To tackle the computational overhead, \cite{cai2010singular} considered singular value thresholding technique. However, it is not scalable to large  problems since it calls singular value decomposition~(SVD) in each iteration. \cite{jaggi2010simple} utilized a sparse semi-definite programming solver to derive a simple yet efficient algorithm. Unfortunately, the memory requirement of their algorithm is proportional to the number of observed entries, making it impractical when the regularized matrix is large and dense (which is the case of LRR). \cite{avron2012efficient} combined stochastic subgradient and incremental SVD to boost efficiency. But for the LRR problem, the type of the loss function does not meet the requirements and thus, it is still not practical to use that algorithm in our case.

Another line in the literature explores a structured formulation of LRR beyond the low-rankness. For example, \cite{wang2013provable} provably showed that combining LRR with SSC can take advantages of both methods. \cite{liu2014recovery} demonstrated that LRR is able to cope with the intrinsic group structure of the data. Very recently, \cite{shen2016learning} argued that the vanilla LRR program does not fully characterize the nature of multiple subspaces, and presented several effective alternatives to LRR.

\subsection{Summary of Contributions}
In this paper, we propose a new algorithm called Online Low-Rank Subspace Clustering (OLRSC), which admits a low computational complexity. In contrast to existing solvers, OLRSC reduces the memory cost of LRR from $\mathcal{O}(n^2)$ to $\mathcal{O}(pd)$~($d < p \ll n$). This nice property makes OLRSC an appealing solution for large-scale subspace clustering problems. Furthermore, we prove that the sequence of solutions produced by OLRSC converges to a stationary point of the expected loss function asymptotically even though only one atom of $\boldsymbol{Y}$ is available at each iteration. In a nutshell, OLRSC resolves {\em all} practical issues of LRR and still promotes global low-rank structure~--~the merit of LRR.

\subsection{Roadmap}
The paper is organized as follows. In Section~\ref{sec:setup}, we reformulate the LRR program~\eqref{eq:lrr} in a way which is amenable for online optimization. Section~\ref{sec:alg} presents the algorithm that incrementally minimizes a surrogate function to the empirical loss. Along with that, we establish a theoretical guarantee in Section~\ref{sec:theory}. The experimental study in Section~\ref{sec:exp} confirms the efficacy and efficiency of our proposed algorithm. Finally, we conclude the work in Section~\ref{sec:conclusion} and the lengthy proof is deferred to the appendix.

\vspace{0.1in}
\noindent{\bf Notation.} \
We use bold lowercase letters, {e.g.} $\boldsymbol{v}$, to denote a column vector. The $\ell_2$ norm and $\ell_1$ norm of a vector $\boldsymbol{v}$ are denoted by $\twonorm{\boldsymbol{v}}$ and $\onenorm{\boldsymbol{v}}$ respectively. Bold capital letters such as $\boldsymbol{M}$ are used to denote a matrix, and its transpose is denoted by $\boldsymbol{M}^\top$. For an invertible matrix $\boldsymbol{M}$, we write its inverse as $\boldsymbol{M}^{-1}$. The capital letter $\boldsymbol{I}_r$ is reserved for identity matrix where $r$ indicates the size. The $j$th column of a matrix $\boldsymbol{M}$ is denoted by $\boldsymbol{m}_j$ if not specified. Three matrix norms  will be used: $\nuclearnorm{\boldsymbol{M}}$ for the nuclear norm, i.e., the sum of the singular values, $\fronorm{\boldsymbol{M}}$ for the Frobenius norm and {$\onenorm{\boldsymbol{M}}$ for the $\ell_1$ norm of a matrix seen as a long vector.} The trace of a square matrix $\boldsymbol{M}$ is denoted by $\tr(\boldsymbol{M})$.
For an integer $n > 0$, we use $[n]$ to denote the integer set $\{1, 2, \cdots, n\}$.

\section{Problem Formulation}\label{sec:setup}
Our goal is to efficiently learn the representation matrix $\boldsymbol{X}$ and the corruption matrix $\boldsymbol{E}$ in an online manner so as to mitigate the issues mentioned in Section~\ref{sec:intro}. The first technique for our purpose is a {\em non-convex reformulation} of the nuclear norm. Assume that the rank of the global optima $\boldsymbol{X}$ in~\eqref{eq:lrr} is at most $d$. Then a standard result in the literature (see, e.g.,~\cite{fazel2001rank}) showed that,
\begin{equation}
\label{eq:nuclear reform}
\nuclearnorm{\boldsymbol{X}} = \min_{\boldsymbol{U}, \boldsymbol{V}, \boldsymbol{X}=\boldsymbol{U} \boldsymbol{V}^\top} \fractwo{1} \( \fronorm{\boldsymbol{U}}^2 + \fronorm{\boldsymbol{V}}^2 \),
\end{equation}
where $\boldsymbol{U} \in \mathbb{R}^{n\times d}$ and $\boldsymbol{V} \in \mathbb{R}^{n\times d}$. The minimum can be attained at, for example, $\boldsymbol{U} = \boldsymbol{U}_0 \boldsymbol{S}_0^{\frac{1}{2}}$ and $\boldsymbol{V} = \boldsymbol{V}_0 \boldsymbol{S}_0^{\frac{1}{2}}$ where $\boldsymbol{X} = \boldsymbol{U}_0 \boldsymbol{S}_0 \boldsymbol{V}_0^\top$ is the singular value decomposition.

In this way, \eqref{eq:lrr} can be written as follows:
\begin{align}
\min_{\boldsymbol{U}, \boldsymbol{V}, \boldsymbol{E}}\ \fractwo{\lambda_1}\fronorm{\boldsymbol{Z}-\boldsymbol{Y}\boldsymbol{U}\boldsymbol{V}^\top-\boldsymbol{E}}^2 + \fractwo{1} \fronorm{\boldsymbol{U}}^2  + \fractwo{1}\fronorm{\boldsymbol{V}}^2  + \lambda_2 \onenorm{\boldsymbol{E}}.
\end{align}
Note that by this reformulation, updating the entries in $\boldsymbol{X}$ amounts to sequentially updating the rows of $\boldsymbol{U}$ and $\boldsymbol{V}$. Also note that this technique is utilized in~\cite{feng2013online} for online RPCA.  Unfortunately, the size of $\boldsymbol{U}$ and $\boldsymbol{V}$ in our problem are both proportional to $n$ and the dictionary $\boldsymbol{Y}$ is partially observed in each iteration, making the algorithm in~\cite{feng2013online} not applicable to LRR. Related to online implementation, another challenge is that, all the rows of $\boldsymbol{U}$ are coupled together at this moment as $\boldsymbol{U}$ is left multiplied by $\boldsymbol{Y}$ in the first term. This makes it difficult to sequentially update the rows of $\boldsymbol{U}$.

For the sake of decoupling the rows of $\boldsymbol{U}$, as part of the crux of our technique, we introduce an {auxiliary variable} $\boldsymbol{D} = \boldsymbol{Y}\boldsymbol{U}$, whose size is $p \times d$ ({i.e.}, independent of the sample size $n$). Interestingly, in this way, we are approximating  the term $\boldsymbol{Z} - \boldsymbol{E}$  with $\boldsymbol{D}\boldsymbol{V}^\top$, which provides an intuition on the role of $\boldsymbol{D}$: Namely, $\boldsymbol{D}$ can be seen as a {\em basis dictionary} of the clean data, with $\boldsymbol{V}$ being the coefficients.

These key observations allow us to derive an equivalent reformulation of LRR~\eqref{eq:lrr}:
\begin{align}
\min_{\boldsymbol{D}, \boldsymbol{U}, \boldsymbol{V}, \boldsymbol{E}}\ \fractwo{\lambda_1}\fronorm{\boldsymbol{Z}-\boldsymbol{Y}\boldsymbol{U}\boldsymbol{V}^\top-\boldsymbol{E}}^2 + \fractwo{1} \( \fronorm{\boldsymbol{U}}^2 + \fronorm{\boldsymbol{V}}^2 \)  + \lambda_2 \onenorm{\boldsymbol{E}},\quad \mathrm{s.t.}\ \ \boldsymbol{D} = \boldsymbol{Y}\boldsymbol{U}.
\end{align}
By penalizing the constraint in the objective, we obtain a {\em regularized} version of LRR on which our algorithm is based
\begin{align}
\label{eq:reg batch}
\min_{\boldsymbol{D}, \boldsymbol{U}, \boldsymbol{V}, \boldsymbol{E}}\ \fractwo{\lambda_1} \fronorm{ \boldsymbol{Z} - \boldsymbol{D}\boldsymbol{V}^\top -\boldsymbol{E} }^2 + \fractwo{1}\( \fronorm{\boldsymbol{U}}^2 + \fronorm{\boldsymbol{V}}^2 \)  + \lambda_2 \onenorm{\boldsymbol{E}} + \fractwo{\lambda_3} \fronorm{\boldsymbol{D} - \boldsymbol{Y}\boldsymbol{U}}^2.
\end{align}
\begin{remark}[Superiority to LRR]
There are two advantages of \eqref{eq:reg batch} compared to~\eqref{eq:lrr}. First, it is amenable for online optimization. Second, it is more informative since it explicitly models the basis of the union of subspaces, hence a better subspace recovery and clustering (see Section~\ref{sec:exp}). This actually meets the core idea of~\cite{liu2014recovery} but they {\em assumed} $Y$ contains true subspaces whereas we {\em learn} the true subspaces.
\end{remark}
\begin{remark}[Parameter]
Note that $\lambda_3$ may be gradually increased until some maximum value is attained so as to enforce the equality constraint. In this way, \eqref{eq:reg batch} attains the same minimum as \eqref{eq:lrr}. Actually, the choice of $\lambda_3$ depends on how much information $\boldsymbol{Y}$ brings for the subspace basis. As we aforementioned, $\boldsymbol{D}$ is the basis dictionary of the clean data and is in turn approximated by (or equal to) $\boldsymbol{Y}\boldsymbol{U}$. This suggests that the range of $\boldsymbol{D}$ is a subset of that of $\boldsymbol{Y}$. As a typical choice of $\boldsymbol{Y} = \boldsymbol{Z}$, if $\boldsymbol{Z}$ is slightly corrupted, we would like to pick a large quantity for $\lambda_3$.
\end{remark}
\begin{remark}[Connection to RPCA]
Due to our explicit modeling of the basis, we unify  LRR and RPCA as follows: for LRR, $\boldsymbol{D} \approx \boldsymbol{Y}\boldsymbol{U}$ (or $\boldsymbol{D} = \boldsymbol{Y}\boldsymbol{U}$ if $\lambda_3$ tends to infinity) while for RPCA, $\boldsymbol{D} = \boldsymbol{U}$. That is, ORPCA~\cite{feng2013online} considers a problem of $\boldsymbol{Y} = \boldsymbol{I}_p$ whose size is independent of $n$, hence can be kept in memory which naturally resolves Issue~\ref{is:memY} and~\ref{is:partialY}. This is why RPCA can be easily implemented in an online fashion while LRR cannot.
\end{remark}
\begin{remark}[Connection to Dictionary Learning]
Generally speaking, LRR~\eqref{eq:lrr} can be seen as a coding algorithm, with the dictionary $\boldsymbol{Y}$ known in advance and $\boldsymbol{X}$ is a desired structured code while other popular algorithms such as dictionary learning~(DL)~\cite{mairal2010online} simultaneously optimizes the dictionary and the sparse code. Interestingly, in view of~\eqref{eq:reg batch}, the link of LRR and DL becomes more clear in the sense that the difference lies in the way how the dictionary is constrained. That is, for LRR we have $\boldsymbol{D} \approx \boldsymbol{Y}\boldsymbol{U}$ and $\boldsymbol{U}$ is further regularized by Frobenius norm whereas for DL, we have $\twonorm{\boldsymbol{d}_i} \leq 1$ for each column of $D$.
\end{remark}

Let $\boldsymbol{z}_i$, $\boldsymbol{y}_i$, $\boldsymbol{e}_i$, $\boldsymbol{u}_i$, and $\boldsymbol{v}_i$ be the $i$th column of matrices $\boldsymbol{Z}$, $\boldsymbol{Y}$, $\boldsymbol{E}$, $\boldsymbol{U}^\top$ and $\boldsymbol{V}^\top$ respectively and define
\begin{align}
\label{eq:tl}
\tilde{\ell}(\boldsymbol{z}, \boldsymbol{D}, \boldsymbol{v}, \boldsymbol{e}) &\stackrel{\text{def}}{=} \fractwo{\lambda_1} \twonorm{ \boldsymbol{z} - \boldsymbol{D} \boldsymbol{v} - \boldsymbol{e} }^2 + \fractwo{1}\twonorm{\boldsymbol{v}}^2 + \lambda_2 \onenorm{\boldsymbol{e}},\\
\label{eq:ell}
\ell(\boldsymbol{z}, \boldsymbol{D}) &= \min_{\boldsymbol{v}, \boldsymbol{e}} \tilde{\ell}(\boldsymbol{z}, \boldsymbol{D}, \boldsymbol{v}, \boldsymbol{e}).
\end{align}
In addition, let
\begin{align}
\label{eq:th}
\th(\boldsymbol{Y}, \boldsymbol{D}, \boldsymbol{U}) &\stackrel{\text{def}}{=} \sum_{i=1}^n \fractwo{1} \twonorm{\boldsymbol{u}_i}^2 + \fractwo{\lambda_3} \fronorm{\boldsymbol{D} - \sum_{i=1}^{n} \boldsymbol{y}_i \boldsymbol{u}_i^\top}^2,\\
\label{eq:h}
h(\boldsymbol{Y}, \boldsymbol{D}) &= \min_{\boldsymbol{U}} \th(\boldsymbol{Y},\boldsymbol{D},\boldsymbol{U}).
\end{align}
Then \eqref{eq:reg batch} can be rewritten as:
\begin{align}\label{eq:tmp}
\min_{\boldsymbol{D}} \min_{\boldsymbol{U}, \boldsymbol{V}, \boldsymbol{E}}\ \sum_{i=1}^{n} \tilde{\ell}(\boldsymbol{z}_i, \boldsymbol{D}, \boldsymbol{v}_i, \boldsymbol{e}_i) + \th(\boldsymbol{Y}, \boldsymbol{D}, \boldsymbol{U}),
\end{align}
which amounts to minimizing the empirical loss function:
\begin{align}
\label{eq:f_n(D)}
f_n(\boldsymbol{D}) \stackrel{\text{def}}{=} \frac{1}{n}\sum_{i=1}^{n} \ell(\boldsymbol{z}_i, \boldsymbol{D}) + \frac{1}{n} h(\boldsymbol{Y}, \boldsymbol{D}).
\end{align}
{In stochastic optimization, we are interested in analyzing the optimality of the obtained solution with respect to the expected loss function. To this end, we first derive the optimal solutions $\boldsymbol{U}^*$, $\boldsymbol{V}^*$ and $\boldsymbol{E}^*$ that minimize~\eqref{eq:tmp} which renders a concrete form of the empirical loss function $f_n(\boldsymbol{D})$, hence we are able to derive the expected loss.

Given $\boldsymbol{D}$, we need to compute the optimal solutions $\boldsymbol{U}^*$, $\boldsymbol{V}^*$ and $\boldsymbol{E}^*$ to evaluate the objective value of $f_n(\boldsymbol{D})$. What is of interest here is that, the optimization procedure of $\boldsymbol{U}$ is totally different from that of $\boldsymbol{V}$ and $\boldsymbol{E}$. According to \eqref{eq:ell}, when $\boldsymbol{D}$ is given, each $\boldsymbol{v}_i^*$ and $\boldsymbol{e}_i^*$ can be solved by only accessing the $i$th sample $\boldsymbol{z}_i$. However, the optimal $\boldsymbol{u}_i^*$ depends on the whole dictionary $\boldsymbol{Y}$ as the second term in $\th(\boldsymbol{Y}, \boldsymbol{D}, \boldsymbol{U})$ couples all the $\boldsymbol{u}_i$'s. Fortunately, it is possible to obtain a closed form solution for $\boldsymbol{U}^*$ which simplifies our analysis. To be more precise, the first order optimality condition for~\eqref{eq:h} gives
\begin{align}
\frac{\partial \th(\boldsymbol{Y}, \boldsymbol{D}, \boldsymbol{U})}{\partial \boldsymbol{U}} =\ \boldsymbol{U} + \lambda_3 (\boldsymbol{Y}^\top \boldsymbol{Y}\boldsymbol{U} - \boldsymbol{Y}^\top \boldsymbol{D}) = 0,
\end{align}
which implies
\begin{align}
\label{eq:U^*}
\boldsymbol{U}^* &= \({\lambda_3}^{-1}\boldsymbol{I}_p + \boldsymbol{Y}^\top \boldsymbol{Y}\)^{-1} \boldsymbol{Y}^\top \boldsymbol{D} \notag\\
&= \lambda_3 \sum_{j=0}^{+\infty} \(-\lambda_3 \boldsymbol{Y}^\top \boldsymbol{Y}\)^j \boldsymbol{Y}^\top \boldsymbol{D} \notag\\
&= \lambda_3 \boldsymbol{Y}^\top \Bigg[ \sum_{j=0}^{+\infty} \(-\lambda_3 \boldsymbol{Y} \boldsymbol{Y}^\top\)^j \Bigg] \boldsymbol{D} \notag\\
&=\boldsymbol{Y}^\top \({\lambda_3}^{-1}\boldsymbol{I}_p + \boldsymbol{Y}\bY^\top\)^{-1} \boldsymbol{D}.
\end{align}
Likewise, another component $\boldsymbol{Y}\boldsymbol{U}^{*\top}$ in~\eqref{eq:th} can be derived as follows:
\begin{align}
\boldsymbol{Y}\boldsymbol{U}^{*\top} = \boldsymbol{D} - \frac{1}{n} \(\frac{1}{n} \boldsymbol{I}_p+\frac{\lambda_3}{n} \boldsymbol{N}_n \)^{-1}\boldsymbol{D},
\end{align}
where we denote
\begin{align}
\label{eq:N_n}
\boldsymbol{N}_n = \sum_{i=1}^{n} \boldsymbol{y}_i \boldsymbol{y}_i^\top.
\end{align}
Recall that $\boldsymbol{u}_i$ is the $i$th column of $\boldsymbol{U}^\top$. So for each $i \in [n]$, we immediately have
\begin{align}
\label{eq:actual u_t}
\boldsymbol{u}_i^* = \boldsymbol{D}^\top \(\frac{1}{\lambda_3}\boldsymbol{I}_p + \sum_{i=1}^{n} \boldsymbol{y}_i \boldsymbol{y}^\top_i\)^{-1} \boldsymbol{y}_i = \frac{1}{n} \boldsymbol{D}^\top \(\frac{1}{\lambda_3 n}\boldsymbol{I}_p + \frac{1}{n} \boldsymbol{N}_n\)^{-1} \boldsymbol{y}_i.
\end{align}
Plugging $\boldsymbol{U}^*$ and $\boldsymbol{Y}\boldsymbol{U}^{*\top}$ back to $\th(\boldsymbol{Y}, \boldsymbol{D}, \boldsymbol{U})$ gives
\begin{align}
h(\boldsymbol{Y}, \boldsymbol{D}) = \frac{1}{n^2} \sum_{i=1}^n \fractwo{1} \twonorm{  \boldsymbol{D}^\top \(\frac{1}{\lambda_3 n}\boldsymbol{I}_p + \frac{1}{n} \boldsymbol{N}_n\)^{-1} \boldsymbol{y}_i }^2  + \frac{\lambda_3}{2 n^2} \fronorm{ \(\frac{1}{n} \boldsymbol{I}_p +\frac{\lambda_3}{n} \boldsymbol{N}_n \)^{-1}\boldsymbol{D} }^2.
\end{align}

Now we derive the expected loss function, which is defined as the limit of the empirical loss function when $n$ tends to infinity. If we assume that all the samples are drawn independently and identically from some (unknown) distribution, we have
\begin{align}
\lim_{n \rightarrow \infty} \frac{1}{n} \sum_{i=1}^{n} \ell(\boldsymbol{z}_i, \boldsymbol{D}) = \mathbb{E}_{\boldsymbol{z}} [\ell(\boldsymbol{z}, \boldsymbol{D})].
\end{align}
If we further assume that the smallest singular value of $\frac{1}{n} \boldsymbol{N}_n$ is bounded away from zero~(which implies $\boldsymbol{N}_n$ is invertible and the spectrum of $\boldsymbol{N}_n^{-1}$ is bounded from the above), we have
\begin{align}
0 \leq \lim_{n \rightarrow \infty} \frac{1}{n} h(\boldsymbol{Y},\boldsymbol{D})  \leq \lim_{n \rightarrow \infty} \frac{1}{n^3} \sum_{i=1}^n \mathrm{C}_0 = 0.
\end{align}
Here $\mathrm{C}_0$ is some absolute constant since $\boldsymbol{D}$ is fixed and $\boldsymbol{y}_i$'s are bounded. Hence, it follows that
\begin{align}
\lim_{n \rightarrow \infty} \frac{1}{n} h(\boldsymbol{Y},\boldsymbol{D}) = 0.
\end{align}
Finally, the expected loss function is given by
\begin{align}
\label{eq:f(D)}
f(\boldsymbol{D}) \stackrel{\text{def}}{=} \lim_{n \rightarrow \infty} f_n(\boldsymbol{D}) = \mathbb{E}_{\boldsymbol{z}}\big[\ell(\boldsymbol{z}, \boldsymbol{D}) \big].
\end{align}


\section{Algorithm}\label{sec:alg}
Our OLRSC algorithm is summarized in Algorithm~\ref{alg:all}. Recall that OLRSC is an online implementation to solve \eqref{eq:f_n(D)}, which is derived from the regularized version of LRR~\eqref{eq:reg batch}. The main idea is optimizing the variables in an alternative manner. That is, at the $t$-th iteration, assume the basis dictionary $\boldsymbol{D}_{t-1}$ is given, we compute the optimal solutions $\{\boldsymbol{v}_t, \boldsymbol{e}_t\}$ by minimizing the objective function $\tilde{\ell}(\boldsymbol{z}_t, {\boldsymbol{D}}_{t-1}, \boldsymbol{v}, \boldsymbol{e})$ over $\boldsymbol{v}$ and $\boldsymbol{e}$. For $\boldsymbol{u}_t$, we need a more carefully designed paradigm since a direct optimization involves loading the full dictionary $\boldsymbol{Y}$ (see~\eqref{eq:actual u_t}). We will elaborate the details later. Subsequently, we update the basis dictionary $\boldsymbol{D}_{t}$ by optimizing a surrogate function to the empirical loss $f_n(\boldsymbol{D})$. In our algorithm, we need to maintain three additional accumulation matrices for which the sizes are independent of $n$.


\begin{algorithm}[t]
\caption{Online Low-Rank Subspace Clustering}
\label{alg:all}
\begin{algorithmic}[1]
    \REQUIRE $\boldsymbol{Z} \in \mathbb{R}^{p\times n}$ (observed samples), $\boldsymbol{Y} \in \mathbb{R}^{p\times n}$, parameters $\lambda_1$, $\lambda_2$ and $\lambda_3$, initial basis $\boldsymbol{D}_0 \in \mathbb{R}^{p\times d}$, zero matrices $\boldsymbol{M}_0 \in \mathbb{R}^{p\times d}$, $\boldsymbol{A}_0 \in \mathbb{R}^{d\times d}$ and $\boldsymbol{B}_0 \in \mathbb{R}^{p\times d}$.
    \ENSURE Optimal basis $\boldsymbol{D}_n$.
    \FOR{$t=1$ to $n$}
      \STATE Access the $t$-th sample $\boldsymbol{z}_t$ and the $t$-th atom $\boldsymbol{y}_t$.
      \STATE Compute the coefficient and noise:
        \begin{align*}
         & \{\boldsymbol{v}_t, \boldsymbol{e}_t\} = \argmin_{\boldsymbol{v}, \boldsymbol{e}} \tilde{\ell}({\boldsymbol{z}_t, \boldsymbol{D}_{t-1}, \boldsymbol{v}, \boldsymbol{e}}),\\
          &\boldsymbol{u}_t = \argmin_{\boldsymbol{u}} \tilde{\ell}_2(\boldsymbol{y}_t, \boldsymbol{D}_{t-1}, \boldsymbol{M}_{t-1}, \boldsymbol{u}).
        \end{align*}\vspace{-0.2in}
      \STATE Update the accumulation matrices:
          \begin{align*}
          \boldsymbol{M}_t \leftarrow \boldsymbol{M}_{t-1} + \boldsymbol{y}_t \boldsymbol{u}_t^\top,\ \ \boldsymbol{A}_t \leftarrow \boldsymbol{A}_{t-1} + \boldsymbol{v}_t\boldsymbol{v}_t^\top,\ \ \boldsymbol{B}_t \leftarrow \boldsymbol{B}_{t-1} + (\boldsymbol{z}_t - \boldsymbol{e}_t)\boldsymbol{v}_t^\top.
          \end{align*}\vspace{-0.2in}
      \STATE Update the basis dictionary:
      	  \begin{align*}
      	 \boldsymbol{D}_t = \argmin_{\boldsymbol{D} } \frac{1}{t} \Bigg[ \frac{1}{2}\tr\(\boldsymbol{D}^\top \boldsymbol{D}(\lambda_1 \boldsymbol{A}_t + \lambda_3 \boldsymbol{I}_d)\)  - \tr\(\boldsymbol{D}^\top (\lambda_1 \boldsymbol{B}_t + \lambda_3 \boldsymbol{M}_t)\) \Bigg].
      	  \end{align*}
    \ENDFOR
\end{algorithmic}
\end{algorithm}

\vspace{0.1in}
\noindent{\bf Solving $\{\boldsymbol{v}_t, \boldsymbol{e}_t\}$.} \ 
We observe that if $\boldsymbol{e}$ is fixed, we can optimize $\boldsymbol{v}$ in closed form:
\begin{align}\label{eq:v}
\boldsymbol{v} = \( \lambda_1^{-1}\boldsymbol{I}_d + \boldsymbol{D}_{t-1}^\top \boldsymbol{D}_{t-1}\)^{-1}\boldsymbol{D}_{t-1}^\top \(\boldsymbol{z}_t - \boldsymbol{e}\).
\end{align}
Conversely, given $\boldsymbol{v}$, the variable $\boldsymbol{e}$ is obtained via soft-thresholding~\cite{donoho1995denoising}:
\begin{align}\label{eq:e}
\boldsymbol{e} = \mathcal{S}_{\lambda_2 / \lambda_1}\(\boldsymbol{z}_t - \boldsymbol{D}_{t-1}\boldsymbol{v}\).
\end{align}
Thus, we utilize block coordinate minimization algorithm to optimize $\boldsymbol{v}$ and $\boldsymbol{e}$.

\vspace{0.1in}
\noindent{\bf Solving $\boldsymbol{u}_t$.} \ 
The closed form solution~\eqref{eq:actual u_t} tells us that it is impossible to derive an accurate estimation of $\boldsymbol{u}_t$ without the entire dictionary $\boldsymbol{Y}$. Thus, we have to ``approximately'' solve it during the online optimization procedure\footnote{Here, ``accurately'' and ``approximately'' mean that when only $\boldsymbol{D}_{t-1}$, $\boldsymbol{z}_t$ and $\boldsymbol{y}_t$ are given, whether we can obtain the same solution $\{\boldsymbol{v}_t, \boldsymbol{e}_t, \boldsymbol{u}_t\}$ as for the batch problem~\eqref{eq:f_n(D)}.}.

Our carefully designed approximate process to solve $\th(\boldsymbol{Y}, \boldsymbol{D}, \boldsymbol{U})$~\eqref{eq:th} is motivated by the coordinate minimization method appealing to $\th(\boldsymbol{Y}, \boldsymbol{D}, \boldsymbol{U})$. As a convention, such method starts with initial guess that $\boldsymbol{u}_i = \boldsymbol{0}$ for all $i \in [n]$ and updates the $\boldsymbol{u}_i$'s in a cyclic order, i.e., $\boldsymbol{u}_1$, $\boldsymbol{u}_2$, $\cdots$, $\boldsymbol{u}_n$, $\boldsymbol{u}_1$, $\cdots$. Let us consider the first pass where we have already updated $\boldsymbol{u}_1$, $\boldsymbol{u}_2$, $\cdots$, $\boldsymbol{u}_{t-1}$ and are to optimize over $\boldsymbol{u}_t$ for some $t > 0$. Note that since the initial values are zero, $\boldsymbol{u}_{t+1} = \boldsymbol{u}_{t+2} = \cdots = \boldsymbol{u}_n = \boldsymbol{0}$. Thereby, the optimal $\boldsymbol{u}_t$ is actually given by minimizing the following function:
\begin{align}
\label{eq:tl_2}
\tilde{\ell}_2(\boldsymbol{y}_t, \boldsymbol{D}, \boldsymbol{M}_{t-1}, \boldsymbol{u})  \stackrel{\text{def}}{=} \fractwo{1} \twonorm{\boldsymbol{u}}^2 + \fractwo{\lambda_3} \fronorm{\boldsymbol{D} - \boldsymbol{M}_{t-1} - \boldsymbol{y}_t \boldsymbol{u}^\top}^2,
\end{align}
where
\begin{align}
\boldsymbol{M}_{t-1} = \sum_{i=1}^{t-1}\boldsymbol{y}_i \boldsymbol{u}_i^\top.
\end{align}
We easily obtain the closed form solution to~\eqref{eq:tl_2} as follows:
\begin{align}
\label{eq:u_t}
\boldsymbol{u}_t = (\twonorm{\boldsymbol{y}_t}^2 + {1}/{\lambda_3})^{-1} (\boldsymbol{D} - \boldsymbol{M}_{t-1})^\top \boldsymbol{y}_t.
\end{align}

Now let us turn to the alternating minimization algorithm, where $\boldsymbol{D}$ is updated iteratively rather than fixed as in~\eqref{eq:u_t}. The above coordinate minimization process can be adjusted in this scenario as we did in Algorithm~\ref{alg:all}. That is, given $\boldsymbol{D}_{t-1}$, after revealing a new atom $\boldsymbol{y}_t$, we compute $\boldsymbol{u}_t$ by minimizing $\tilde{\ell}_2(\boldsymbol{y}_t, \boldsymbol{D}_{t-1}, \boldsymbol{M}_{t-1}, \boldsymbol{u})$, followed by updating $\boldsymbol{D}_t$. In this way, when the algorithm terminates, we in essence run a one-pass update on $\boldsymbol{u}_t$'s with a simultaneous computation of new basis dictionary.



\vspace{0.1in}
\noindent{\bf Solving $\boldsymbol{D}_t$.} \ 
As soon as the past filtration $\{ \boldsymbol{v}_i, \boldsymbol{e}_i, \boldsymbol{u}_i \}_{i=1}^t$ are available, we can compute a new iterate $\boldsymbol{D}_t$ by optimizing the surrogate function
\begin{align}
\label{eq:g_t(D)}
g_t(\boldsymbol{D}) \stackrel{\text{def}}{=} \frac{1}{t} \(\sum_{i=1}^{t} \tilde{\ell}(\boldsymbol{z}_i, \boldsymbol{D}, \boldsymbol{v}_i, \boldsymbol{e}_i) + \sum_{i=1}^{t} \fractwo{1} \twonorm{\boldsymbol{u}_i}^2  + \fractwo{\lambda_3} \fronorm{\boldsymbol{D} - \boldsymbol{M}_t}^2 \).
\end{align}
Expanding the first term, we find that $\boldsymbol{D}_t$ is given by
\begin{align}
\boldsymbol{D}_t =&\ \argmin_{\boldsymbol{D} } \frac{1}{t} \Bigg[ \frac{1}{2}\tr\(\boldsymbol{D}^\top \boldsymbol{D}(\lambda_1 \boldsymbol{A}_t + \lambda_3 \boldsymbol{I}_d)\)  - \tr\(\boldsymbol{D}^\top (\lambda_1 \boldsymbol{B}_t + \lambda_3 \boldsymbol{M}_t)\) \Bigg]\notag\\
=&\  (\lambda_1 \boldsymbol{B}_t + \lambda_3 \boldsymbol{M}_t)  (\lambda_1 \boldsymbol{A}_t + \lambda_3 \boldsymbol{I}_d)^{-1},
\end{align}
where $\boldsymbol{A}_t = \sum_{i=1}^{t} \boldsymbol{v}_i \boldsymbol{v}_i^\top$ and $\boldsymbol{B}_t = \sum_{i=1}^{t} (\boldsymbol{z}_i - \boldsymbol{e}_i)\boldsymbol{v}_i^\top$. We point out that the size of $\boldsymbol{A}_t$ is $d\times d$ and that of $\boldsymbol{B}_t$ is $p\times d$, i.e., independent of sample size. In practice, as~\cite{mairal2010online} suggested, one may apply a block coordinate descent approach to minimize over $\boldsymbol{D}$. Compared to the closed form solution given above, such algorithm usually converges very fast after revealing sufficient number of samples. In fact, we observe that a one-pass update on the columns of $\boldsymbol{D}$ suffices to ensure a favorable performance. See Algorithm~\ref{alg:D}.

\begin{algorithm}[h]
\caption{Solving $\boldsymbol{D}$}
\label{alg:D}
\begin{algorithmic}[1]
    \REQUIRE $\boldsymbol{D} \in \mathbb{R}^{p\times d}$ in the previous iteration, accumulation matrix $\boldsymbol{M}$, $\boldsymbol{A}$ and $\boldsymbol{B}$, parameters $\lambda_1$ and $\lambda_3$.
    \ENSURE Optimal $\boldsymbol{D}$ (updated).
    \STATE Denote $\widehat{\boldsymbol{A}} = \lambda_1 \boldsymbol{A} + \lambda_3 \boldsymbol{I}$ and $\widehat{\boldsymbol{B}} = \lambda_1 \boldsymbol{B} + \lambda_3 \boldsymbol{M} $.
    \REPEAT
      \FOR{$j = 1$ to $d$}
      \STATE Update the $j$th column of $\boldsymbol{D}$:
      \begin{align*}
        \boldsymbol{d}_j \leftarrow \boldsymbol{d}_j - \frac{1}{\widehat{A}_{jj}} \( \boldsymbol{D} \widehat{\boldsymbol{a}}_j - \widehat{\boldsymbol{b}}_j \)
      \end{align*}
      \ENDFOR
    \UNTIL{convergence}
\end{algorithmic}
\end{algorithm}


\vspace{0.1in}
\noindent{\bf Memory Cost.} \
It is remarkable that the  memory cost of Algorithm~\ref{alg:all} is $\mathcal{O}(pd)$. To see this, note that when solving $\boldsymbol{v}_t$ and $\boldsymbol{e}_t$, we load the auxiliary variable $\boldsymbol{D}_t$ and a sample $\boldsymbol{z}_t$ into the memory, which costs $\mathcal{O}(pd)$. To compute the optimal $\boldsymbol{u}_t$'s, we need to access $\boldsymbol{D}_t$ and $\boldsymbol{M}_t \in \mathbb{R}^{p\times d}$. Although we aim to minimize \eqref{eq:g_t(D)}, which seems to require all the past information, we actually only need to record $\boldsymbol{A}_t$, $\boldsymbol{B}_t$ and $\boldsymbol{M}_t$, whose sizes are at most $\mathcal{O}(pd)$ (since $d < p$).

\vspace{0.1in}
\noindent{\bf Computational Efficiency.} \
In addition to  memory efficiency,  we further clarify that the the computation in each iteration is cheap. To compute $\{\boldsymbol{v}_t, \boldsymbol{e}_t\}$, one may utilize the block coordinate method in~\cite{richtarik2014iteration} which enjoys linear convergence due to strong convexity. One may also apply the stochastic variance reduced algorithms which also ensure a geometric rate of convergence~\cite{xiao2014proximal,saga}. The $\boldsymbol{u}_t$ is given by simple matrix-vector multiplication, which costs $\mathcal{O}(pd)$. It is easy to see the update on the accumulation matrices is $\mathcal{O}(pd)$ and that of $\boldsymbol{D}_t$ is $\mathcal{O}(pd^2)$.

\vspace{0.1in}
\noindent{\bf A Fully Online Scheme.} \
Now we have provided a way to (approximately) optimize the LRR problem~\eqref{eq:lrr} in online fashion. Usually, researchers in the literature will take an optional post-processing step to refine the segmentation accuracy, for example, applying spectral clustering on the representation matrix $\boldsymbol{X}$. In this case, one has to collect all the $\boldsymbol{u}_i$'s and $\boldsymbol{v}_i$'s to compute $\boldsymbol{X} = \boldsymbol{U} \boldsymbol{V}^\top$ which again increases the memory cost to $\mathcal{O}(n^2)$. Here, we suggest an alternative scheme which admits $\mathcal{O}(kd)$ memory usage where $k$ is the number of subspaces. The idea is utilizing the well-known $k$-means on $\boldsymbol{v}_i$'s. There are two notable advantages compared to the spectral clustering. First, updating the $k$-means model can be implemented in online manner and the computation is $\mathcal{O}(kd)$ for each iteration. Second, we observe that $\boldsymbol{v}_i$ is actually a robust feature for the $i$th sample. Combining the online $k$-means with Algorithm~\ref{alg:all}, we obtain a fully and efficient online subspace clustering scheme where the memory cost is $\mathcal{O}(pd)$. For the reader's convenience, we summarize this pipeline in Algorithm~\ref{alg:fully}.

\begin{algorithm}[h]
\caption{Fully Online Pipeline for Low-Rank Subspace Clustering}
\label{alg:fully}
\begin{algorithmic}[1]
\REQUIRE $\boldsymbol{Z} \in \mathbb{R}^{p\times n}$ (observed samples), $\boldsymbol{Y} \in \mathbb{R}^{p\times n}$, parameters $\lambda_1$, $\lambda_2$ and $\lambda_3$, initial basis $\boldsymbol{D}_0 \in \mathbb{R}^{p\times d}$, zero matrices $\boldsymbol{M}_0 \in \mathbb{R}^{p\times d}$, $\boldsymbol{A}_0 \in \mathbb{R}^{d\times d}$ and $\boldsymbol{B}_0 \in \mathbb{R}^{p\times d}$, number of clusters $k$, initial centroids $\boldsymbol{C} \in \mathbb{R}^{d\times k}$.
\ENSURE Optimal basis $\boldsymbol{D}_n$, cluster centroids $\boldsymbol{C}$, cluster assignments $\{ o_1, o_2, \cdots, o_n \}$.
\STATE Initialize $r_1 = r_2 = \cdots = r_k = 0$.
\FOR{$t=1$ to $n$}
  \STATE Access the $t$-th sample $\boldsymbol{z}_t$ and the $t$-th atom $\boldsymbol{y}_t$.
  \STATE Compute $\{ \boldsymbol{v}_t, \boldsymbol{e}_t, \boldsymbol{u}_t, \boldsymbol{D}_t \}$ by Algorithm~\ref{alg:all}.
\STATE Compute $o_t = \argmin_{1 \leq j \leq k} \twonorm{\boldsymbol{v}_t - \boldsymbol{c}_j}$.
\STATE Update the $o_t$-th center:
\begin{align*}
&r_{o_t} \leftarrow r_{o_t} + 1,\\
&\boldsymbol{c}_{o_t} \leftarrow \frac{r_{o_t} - 1}{r_{o_t}} \boldsymbol{c}_{o_t} + \frac{1}{r_{o_t}} \boldsymbol{v}_t.
\end{align*}
\ENDFOR
\end{algorithmic}
\end{algorithm}

\vspace{0.1in}
\noindent{\bf An Accurate Online Implementation.} \
Our strategy for solving $\boldsymbol{u}_t$ is based on an approximate routine which resolves Issue~\ref{is:partialY} as well as has a low complexity. Yet, to tackle Issue~\ref{is:partialY}, another potential way is to avoid the variable $\boldsymbol{u}_t$\footnote{We would like to thank the anonymous NIPS 2015 Reviewer for pointing out this potential solution to the online algorithm. Here we explain why this alternative  can be computationally expensive.}. Recall that we derive the optimal solution $\boldsymbol{U}^*$ (provided that $\boldsymbol{D}$ is given) to~\eqref{eq:reg batch} is given by~\eqref{eq:U^*}. Plugging it back to~\eqref{eq:reg batch}, we obtain
\begin{align*}
& \fronorm{\boldsymbol{U}^*}^2 = \tr\( \boldsymbol{D} \boldsymbol{D}^\top \( \boldsymbol{Q}_n - {\lambda_3}^{-1} \boldsymbol{Q}^2_n \) \),\\
&\fronorm{\boldsymbol{D} - \boldsymbol{Y}\boldsymbol{U}^*}^2 = \fronorm{\boldsymbol{D} - {\lambda_3}^{-1} \boldsymbol{Q}_n\boldsymbol{D}}^2,
\end{align*}
where
\begin{align*}
\boldsymbol{Q}_n = \( {\lambda_3}^{-1} \boldsymbol{I}_p + \boldsymbol{N}_n \)^{-1}.
\end{align*}
Here, $\boldsymbol{N}_n$ was given in~\eqref{eq:N_n}. Note that the size of $\boldsymbol{Q}_n$ is $p\times p$. Hence, if we incrementally compute the accumulation matrix $\boldsymbol{N}_t = \sum_{i=1}^{t} \boldsymbol{y}_i \boldsymbol{y}_i^\top$, we can update the variable $\boldsymbol{D}$ in an online fashion. Namely, at $t$-th iteration, we re-define the surrogate function as follows:
\begin{align*}
g_t(\boldsymbol{D}) \stackrel{\text{def}}{=} \frac{1}{t} \Bigg[ \sum_{i=1}^{t} \tilde{\ell}(\boldsymbol{z}_i, \boldsymbol{D}, \boldsymbol{v}_i, \boldsymbol{e}_i) + \fractwo{\lambda_3}  \fronorm{\boldsymbol{D} -\frac{1}{\lambda_3} \boldsymbol{Q}_t\boldsymbol{D}}^2 + \fractwo{1} \tr\( \boldsymbol{D} \boldsymbol{D}^\top \( \boldsymbol{Q}_t - \frac{1}{\lambda_3} \boldsymbol{Q}^2_t \) \) \Bigg].
\end{align*}
Again, by noting the fact that $\tilde{\ell}(\boldsymbol{z}_i, \boldsymbol{D}, \boldsymbol{v}_i, \boldsymbol{e}_i)$ only involves recording $\boldsymbol{A}_t$ and $\boldsymbol{B}_t$, we show that the memory cost is independent of sample size.

While promising since the above procedure avoids the approximate computation, the main shortcoming is computing the inverse of a $p\times p$ matrix in each iteration, hence not efficient. Moreover, as we will show in Theorem~\ref{thm:stationary}, although the $\boldsymbol{u}_t$'s are approximate solutions, we are still guaranteed the convergence of $\boldsymbol{D}_t$.


\section{Theoretical Analysis}\label{sec:theory}
We make three assumptions underlying our analysis.
\begin{assumption}\label{as:z}
The observed data are generated i.i.d.\ from some distribution and there exist constants $\alpha_0$ and $\alpha_1$, such that the conditions $0 < \alpha_0 \leq \twonorm{\boldsymbol{z}} \leq \alpha_1$ and $ \alpha_0 \leq \twonorm{\boldsymbol{y}} \leq \alpha_1$ hold almost surely.
\end{assumption}
\begin{assumption}\label{as:N_t}
The smallest singular value of the matrix $\frac{1}{t}\boldsymbol{N}_t$ is lower bounded away from zero.
\end{assumption}
\begin{assumption}\label{as:g_t(D)}
The surrogate functions $g_t(\boldsymbol{D})$ are strongly convex for all $t \geq 0$.
\end{assumption}

Based on these assumptions, we establish the main theoretical result, justifying the validity of {Algorithm~\ref{alg:all}}.
\begin{theorem}
\label{thm:stationary}
Assume~\ref{as:z},~\ref{as:N_t} and~\ref{as:g_t(D)}. Let $\{\boldsymbol{D}_t\}_{t=1}^{\infty}$ be the sequence of optimal bases produced by Algorithm~\ref{alg:all}. Then, the sequence converges to a stationary point of the expected loss function $f(\boldsymbol{D})$ when $t$ goes to infinity.
\end{theorem}
Note that since the reformulation of the nuclear norm~\eqref{eq:nuclear reform} is non-convex, we can only guarantee that the solution is a stationary point in general~\cite{bertsekas1999nonlinear}. We also remark that OLRSC asymptotically fulfills the first order optimality condition of~\eqref{eq:lrr}. To see this, we follow the proof technique of Prop.3 in~\cite{mardani2015subspace} and let $\boldsymbol{X} = \boldsymbol{U} \boldsymbol{V}^\top$, $W_1 = \boldsymbol{U}\bU^\top$, $W_2 = \boldsymbol{V}\bV^\top$, $\boldsymbol{M}_1 = \boldsymbol{M}_3 = 0.5\boldsymbol{I}$, $\boldsymbol{M}_2 = \boldsymbol{M}_4 = 0.5\lambda_1 \boldsymbol{Y}^\top(\boldsymbol{Y}\boldsymbol{X} + \boldsymbol{E} - \boldsymbol{Z})$. Due to our uniform bound~(Prop.~\ref{prop:bound:veABMDu}), we justify the optimality condition. See the details in~\cite{mardani2015subspace}.

More interestingly, as we mentioned in Section~\ref{sec:alg}, the solution~\eqref{eq:u_t} is not accurate in the sense that it is not equal to that of~\eqref{eq:actual u_t} given $\boldsymbol{D}$. Yet, our theorem asserts that this will not deviate $\{\boldsymbol{D}_t\}_{t\geq 0}$ away from the stationary point. The intuition underlying such amazing phenomenon is that the expected loss function~\eqref{eq:f(D)} is only determined by $\ell(\boldsymbol{z}, \boldsymbol{D})$ which does not involve $\boldsymbol{u}_t$. What is of matter for $\boldsymbol{u}_t$ and $\boldsymbol{M}_t$ is their uniform boundedness and concentration to establish the convergence. Thanks to the carefully chosen function $\tilde{\ell}(\boldsymbol{z}, \boldsymbol{D}, \boldsymbol{M}, \boldsymbol{u})$ and the surrogate function $g_t(\boldsymbol{D})$, we are able to prove the desired property by mathematical induction which is a crucial step in our proof.

In particular, we have the following lemma that facilitates our analysis:
\begin{lemma}
Assume~\ref{as:z} and~\ref{as:N_t} and~\ref{as:g_t(D)}. Let $\{\boldsymbol{M}_t\}_{t\geq 0}$ be the sequence of the matrices produced by Algorithm~\ref{alg:all}. Then, there exists some universal constant $\mathrm{C}_0$, such that for all $t \geq 0$, $\fronorm{M_t} \leq \mathrm{C}_0$.
\end{lemma}
Due to the above lemma, the solution $\boldsymbol{D}_t$ is essentially determined by $\frac{1}{t}\boldsymbol{A}_t$ and $\frac{1}{t}\boldsymbol{B}_t$ when $t$ is a very large quantity since $\frac{1}{t}\boldsymbol{M}_t \rightarrow 0$. We also have a non-asymptotic rate for the numerical convergence of $\boldsymbol{D}_t$ as $\twonorm{\boldsymbol{D}_t - \boldsymbol{D}_{t-1}} = \mathcal{O}(1/t)$. See Appendix~\ref{supp:sec:proof} for more details and a full proof.

\section{Experiments}\label{sec:exp}

Before presenting the empirical results, we first introduce the universal settings used throughout the section.

\vspace{0.1in}
\noindent{\bf Algorithms.} \ 
For the subspace recovery task, we compare our algorithm with state-of-the-art solvers including ORPCA~\cite{feng2013online}, LRR~\cite{liu2013robust} and PCP~\cite{candes2011robust}. For the subspace clustering task, we choose ORPCA, LRR and SSC~\cite{elhamifar2009sparse} as the competitive baselines. Recently, \cite{liu2014recovery} improved the vanilla LRR by utilizing some low-rank matrix for $\boldsymbol{Y}$. We denote this variant of LRR by LRR2 and accordingly, our algorithm equipped with such $\boldsymbol{Y}$ is denoted as OLRSC2.

\vspace{0.1in}
\noindent{\bf Evaluation Metric.} \ 
We evaluate the fitness of the recovered subspaces $\boldsymbol{D}$ (with each column being normalized) and the ground truth $L$ by the Expressed Variance (EV)~\cite{xu2010principal}:
\begin{align}
\text{EV}(\boldsymbol{D}, \boldsymbol{L}) \stackrel{\text{def}}{=} \frac{\tr(\boldsymbol{D}\bD^\top \boldsymbol{L}\bL^\top)}{\tr(\boldsymbol{L}\bL^\top)}.
\end{align}
The value of EV scales between 0 and 1, and a higher value means better recovery.

The performance of subspace clustering is measured by clustering accuracy, which also ranges in the interval $[0, 1]$, and a higher value indicates a more accurate clustering.

\vspace{0.1in}
\noindent{\bf Parameters.} \ 
We set $\lambda_1 = 1$, $\lambda_2 = {1}/{\sqrt{p}}$ and $\lambda_3 = \sqrt{t/p}$, where $t$ is the iteration counter. These settings are actually used in ORPCA. We follow the default parameter setting for the baselines.


\subsection{Subspace Recovery}\label{subsec:recovery}

\noindent{\bf Simulation Data.} \
We use 4 disjoint subspaces $\{\mathcal{S}_k \}_{k=1}^4 \subset \mathbb{R}^p$, whose bases are denoted by $\{ L_k \}_{k=1}^4 \in \mathbb{R}^{p \times d_k}$. The clean data matrix $\bar{\boldsymbol{Z}}_k \in \mathcal{S}_k$ is then produced by $\bar{\boldsymbol{Z}}_k = L_k R_k^\top$, where $R_k \in \mathbb{R}^{n_k \times d_k}$. The entries of $L_k$'s and $R_k$'s are sampled i.i.d.\  from the normal distribution. Finally, the observed data matrix $\boldsymbol{Z}$ is generated by $\boldsymbol{Z} = \bar{\boldsymbol{Z}} + \boldsymbol{E}$, where $\bar{\boldsymbol{Z}}$ is the column-wise concatenation of $\bar{\boldsymbol{Z}}_k$'s followed by a random permutation, $\boldsymbol{E}$ is the sparse corruption whose $\rho$ fraction entries are non-zero and follow an i.i.d.\ uniform distribution over $[-2, 2]$. We independently  conduct each experiment 10 times and report the averaged results.

\vspace{0.1in}
\noindent{\bf Robustness.} \
We illustrate by simulation results that OLRSC can effectively recover the underlying subspaces, confirming that $\boldsymbol{D}_t$ converges to the union of subspaces. For the two online algorithms OLRSC and ORPCA, We compute the EV after revealing all the samples. We examine the performance under different intrinsic dimension $d_k$'s and corruption $\rho$. To be more detailed, the $d_k$'s are varied from $0.01p$ to $0.1p$ with a step size $0.01p$, and the $\rho$ is from 0 to 0.5, with a step size 0.05.

\begin{figure}[h!]
\centering
\includegraphics[width=0.24\linewidth]{diff_n_4000_olrsc.eps}
\includegraphics[width=0.24\linewidth]{diff_n_4000_orpca.eps}
\includegraphics[width=0.24\linewidth]{diff_n_4000_lrr.eps}
\includegraphics[width=0.24\linewidth]{diff_n_4000_pcp.eps}
\caption{{\bf Subspace recovery under different intrinsic dimensions and corruptions.} Brighter is better. We set $p = 100$, $n_k = 1000$ and $d = 4d_k$. LRR and PCP are batch methods. OLRSC consistently outperforms ORPCA and even improves the performance of LRR. Compared to PCP, OLRSC is competitive in most cases and degrades a little for highly corrupted data, possibly due to the number of samples is not sufficient for its convergence.
}
\label{fig:diff_rank_rho}
\end{figure}

The results are presented in Figure~\ref{fig:diff_rank_rho}. The most intriguing observation is that OLRSC as an online algorithm outperforms its batch counterpart LRR! Such improvement may come from the explicit modeling for the basis, which makes OLRSC more informative than LRR. Interestingly,~\cite{guo2014online} also observed that in some situations, an online algorithm can outperform the batch counterpart. To fully understand the rationale behind this phenomenon is an important direction for future research. Notably, OLRSC consistently beats ORPCA (an online version of PCP), which may be the consequence of the fact that OLRSC takes into account that the data are produced by a union of small subspaces. While PCP works well for almost all scenarios, OLRSC degrades a little when addressing difficult cases (high rank and corruption). This is not surprising since Theorem~\ref{thm:stationary} is based on asymptotic analysis and hence, we expect that OLRSC will converge to the true subspace after acquiring more samples.

\vspace{0.1in}
\noindent{\bf Convergence Rate.} \
Now we test on a large dataset to show that our algorithm usually converges to the true subspace faster than ORPCA. We plot the EV curve against the number of samples in Figure~\ref{fig:ev_samples}. Firstly, when equipped with a proper matrix $\boldsymbol{Y}$, OLRSC2 and LRR2 can always produce an exact recovery of the subspace as PCP does. When using the dataset itself for $\boldsymbol{Y}$, OLRSC still converges to a favorable point after revealing all the samples. Compared to ORPCA, OLRSC is more robust and converges much faster for hard cases (see, e.g., $\rho = 0.5$). Again, we note that in such hard cases, OLRSC outperforms LRR, which agrees with the observation in Figure~\ref{fig:diff_rank_rho}.


\begin{figure}[h!]
\centering
{\includegraphics[width=0.24\linewidth]{ev_samples_rho_0.01.eps}}
{\includegraphics[width=0.24\linewidth]{ev_samples_rho_0.3.eps}}
{\includegraphics[width=0.24\linewidth]{ev_samples_rho_0.5.eps}}
{\includegraphics[width=0.24\linewidth]{ev_samples_time.eps}}
\caption{{\bf Convergence rate and time complexity.} A higher EV means better subspace recovery. We set $p = 1000$, $n_k = 5000$, $d_k = 25$ and $d = 100$. OLRSC always converges to or outperforms the batch counterpart LRR. For hard cases, OLRSC converges much faster than ORPCA. Both PCP and LRR2 achieve the best EV value. When equipped with the same dictionary as LRR2, OLRSC2 also well handles the highly corrupted data ($\rho=0.5$). Our methods are more efficient than the competitors but PCP when $\rho$ is small, possibly because PCP utilizes a highly optimized C++ toolkit while ours are written in Matlab.
}
\label{fig:ev_samples}
\end{figure}

\vspace{0.1in}
\noindent{\bf Computational Efficiency.} We also illustrate the time complexity of the algorithms in the last panel of Figure~\ref{fig:ev_samples}. In short, our algorithms (OLRSC and OLRSC2) admit the lowest computational complexity for all cases. One may argue that PCP spends slightly less time than ours for a small $\rho$ (0.01 and 0.1). However, we remark here that PCP utilizes a highly optimized C++ toolkit to boost computation while our algorithms are {fully} written in Matlab. We believe that ours will work more efficiently if properly optimized by, {e.g.,} the blas routine. Another important message conveyed by the figure is that, OLRSC is always being orders of magnitude computationally more efficient than the batch method LRR, as well as producing comparable or even better solution.


\subsection{Subspace Clustering}\label{subsec:clustering}
\noindent{\bf Datasets.} \
We examine the performance for subspace clustering on 5 realistic databases shown in Table~\ref{tb:dataset}, which can be downloaded from the LibSVM website. For MNIST, We randomly select 20000 samples to form MNIST-20K since we find it time consuming to run the batch methods on the entire database.

\begin{table}[h]
	\centering
	\caption{\bf Datasets for subspace clustering.}
	\begin{tabular}{lccc}
		\toprule
		& \#classes & \#samples & \#features \\
		\midrule
		Mushrooms & 2         & 8124      & 112        \\
		DNA       & 3         & 3186      & 180        \\
		Protein   & 3         & 24,387    & 357        \\
		USPS      & 10        & 9298      & 256        \\
		MNIST-20K & 10        & 20,000    & 784       \\
		\bottomrule
	\end{tabular}
	\label{tb:dataset}
\end{table}

\vspace{0.1in}
\noindent{\bf Standard Clustering Pipeline.} \
In order to focus on the solution quality of different algorithms, we follow the standard pipeline which feeds $\boldsymbol{X}$ to a spectral clustering algorithm~\cite{ng2002spectral}. To this end, we collect all the $\boldsymbol{u}$'s and $\boldsymbol{v}$'s produced by OLRSC to form the representation matrix $\boldsymbol{X}=\boldsymbol{U}\boldsymbol{V}^\top$. For ORPCA, we use $\boldsymbol{R}_0\boldsymbol{R}_0^\top$ as the similarity matrix~\cite{liu2013robust}, where $\boldsymbol{R}_0$ is the row space of $\boldsymbol{Z}_0=\boldsymbol{L}_0\boldsymbol{\Sigma}_0\boldsymbol{R}_0^\top$ and $\boldsymbol{Z}_0$ is the clean matrix recovered by ORPCA. We run our algorithm and ORPCA with 2 epochs so as to apply backward correction on the coefficients ($\boldsymbol{U}$ and $\boldsymbol{V}$ in ours and $\boldsymbol{R}_0$ in ORPCA).

\vspace{0.1in}
\noindent{\bf Fully Online Pipeline.} \
As we discussed in Section~\ref{sec:alg}, the (optional) spectral clustering procedure needs the similarity matrix $\boldsymbol{X}$, making the memory proportional to $n^2$. To tackle this issue, we proposed a fully online scheme where the key idea is performing $k$-means on $\boldsymbol{V}$. Here, we examine the efficacy of this variant, which is called OLRSC-F.


\begin{table}[h]
\caption{{\bf Clustering accuracy (\%) and computational time (seconds).} For each dataset, the first row indicates the accuracy and the second row the running time. For all the large-scale datasets, OLRSC (or OLRSC-F) has the highest clustering accuracy. Regarding the running time, our method spends comparable time as ORPCA~(the fastest solver) does while dramatically improves the accuracy. Although SSC is slightly better than SSC on Protein, it consumes one hour while OLRSC takes 25 seconds.
}
\centering
\begin{tabular}{lccccc}
\toprule
       & OLRSC &OLRSC-F  & ORPCA & LRR   & SSC \\
\midrule
Mush- & {85.09} & {\bf 89.36}  & 65.26 & 58.44  & 54.16 \\
rooms  & 8.78 & 8.78  & 8.30 & 46.82  & 32 min\\
\midrule
\multirow{2}{*}{DNA}  & {67.11} & {\bf 83.08}  & 53.11 & 44.01  & 52.23\\
 & 2.58 & 2.58 & 2.09 & 23.67 & 3 min\\
\midrule
\multirow{2}{*}{Protein}  & { 43.30} & 43.94  & 40.22 & 40.31  & {\bf 44.27}\\
 & 24.66 & 24.66 & 22.90 & 921.58 & 65 min\\
\midrule
\multirow{2}{*}{USPS}  & {65.95} & {\bf 70.29}  & 55.70 & 52.98  & 47.58\\
 & 33.93 & 33.93  & 27.01 & 257.25 & 50 min\\
\midrule
MNIST-  & {\bf 57.74} & 55.50  & 54.10 & 55.23 &  43.91\\
20K & 129 & 129  & 121 & 32 min &  7 hours\\
\bottomrule
\end{tabular}
\label{tb:sc}
\end{table}

The results are recorded in Table~\ref{tb:sc}, where the time cost of spectral clustering or $k$-means is not included so we can focus on comparing the efficiency of the algorithms themselves. Also note that we use the dataset itself as the dictionary $\boldsymbol{Y}$ because we find that an alternative choice of $\boldsymbol{Y}$ does not help much on this task. For OLRSC and ORPCA, they require an estimation on the true rank. Here, we use $5k$ as such estimation where $k$ is the number of classes of a dataset. Our algorithm significantly outperforms the two state-of-the-art methods LRR and SSC both for accuracy and efficiency. One may argue that SSC is slightly better than OLRSC on Protein. Yet, it spends 1 hour while OLRSC only costs 25 seconds. Hence, SSC is not practical. Compared to ORPCA, OLRSC always identifies more correct samples as well as consumes comparable running time. For example, on the USPS dataset, OLRSC achieves the accuracy of 65.95\% while that of ORPCA is 55.7\%. Regarding the running time, OLRSC uses only 7 seconds more than ORPCA~--~same order of computational complexity, which agrees with the qualitative analysis in Section~\ref{sec:alg}.



\begin{table}[h]
\caption{{\bf Time cost (seconds) of spectral clustering and $k$-means.}}
\centering
\begin{tabular}{lccccc}
\toprule
       & Mushrooms & DNA & Protein & USPS & MNIST-20K\\
\midrule
Spectral & 295 & 18 & 7567 & 482 & 4402\\
$k$-means & 2 & 6 & 5 & 19 & 91\\
\bottomrule
\end{tabular}
\label{tab:fully_online}
\end{table}
More interestingly, it shows that the $k$-means alternative (OLRSC-F) usually outperforms the spectral clustering pipeline. This suggests that perhaps for {\em robust} subspace clustering formulations, the simple $k$-means paradigm suffices to guarantee an appealing result. On the other hand, we report the running time of spectral clustering and $k$-means in Table~\ref{tab:fully_online}. As expected, since spectral clustering computes SVD for an $n$-by-$n$ similarity matrix, it is quite slow. In fact, it sometimes dominates the running time of the whole pipeline. In contrast, $k$-means is extremely fast and scalable, as it can be implemented in online fashion.



\subsection{Influence of $d$}

A key ingredient of our formulation is a factorization on the nuclear norm regularized matrix, which requires an estimation on the rank of the $\boldsymbol{X}$ (see \eqref{eq:nuclear reform}). Here we examine the influence of the selection of $d$ (which plays as an upper bound of the true rank). We report both EV and clustering accuracy for different $d$ under a range of corruptions. The simulation data are generated as in Section~\ref{subsec:recovery} and we set $p=200$, $n_k=1000$ and $d_k = 10$. Since the four subspaces are disjoint, the true rank is 40.

\begin{figure}[h]
\centering
	{\includegraphics[width=0.3\linewidth]{diff_d_EV.eps}}
	{\includegraphics[width=0.3\linewidth]{diff_d_ClusteringAccuracy.eps}}
	\caption{{\bf Examine the influence of ${d}$.} We experiment on $d = \{ 2, 20, 40, 60, 80, 100, 120, 140, 160, 180 \}$.  The true rank is 40. }
	\label{fig:diff_d}
\end{figure}

From Figure~\ref{fig:diff_d}, we observe that our algorithm cannot recover the true subspace if $d$ is smaller than the true rank. On the other hand, when $d$ is sufficiently large (at least larger than the true rank), our algorithm can perfectly estimate the subspace. This agrees with the results in~\cite{burer2005local} which says as long as $d$ is large enough, any local minima is global optima. We also illustrate the influence of $d$ on subspace clustering. Generally speaking, OLRSC can consistently identify the cluster of the data points if $d$ is sufficiently large. Interestingly, different from the subspace recovery task, here the requirement for $d$ seems to be slightly relaxed. In particular, we notice that if we pick $d$ as 20 (smaller than the true rank), OLRSC still performs well. Such relaxed requirement of $d$ may benefit from the fact that the spectral clustering step can correct some wrong points as suggested by~\cite{soltanolkotabi2014robust}.


\section{Conclusion}\label{sec:conclusion}
In this paper, we have proposed an online algorithm called OLRSC for subspace clustering, which dramatically reduces the memory cost of LRR from $\mathcal{O}(n^2)$ to $\mathcal{O}(pd)$~--~ orders of magnitudes more memory efficient. One of the key techniques in this work is an explicit basis modeling, which essentially renders the model more informative than LRR. Another important component is a non-convex reformulation of the nuclear norm. Combining these techniques allows OLRSC to simultaneously recover the union of the subspaces, identify the possible corruptions and perform subspace clustering. We have also established the theoretical guarantee that solutions produced by our algorithm  converge to a stationary point of the expected loss function. Moreover, we have analyzed the time complexity and empirically demonstrated that our algorithm is computationally very efficient compared to competing baselines. Our extensive experimental study on synthetic and realistic datasets  also illustrates the robustness of OLRSC. In a nutshell, OLRSC is an appealing algorithm  in all three worlds: memory cost, computation and robustness.

\clearpage



