Online Robust Subspace Tracking from Partial Information

Introduction

Low-rank subspaces have long been a powerful tool in data modeling and analysis. Applications in communications , source localization and target tracking in radar and sonar , and medical imaging all leverage subspace models in order to recover the signal of interest and reject noise. In these classical signal processing problems, a handful of high-quality sensors are co-located such that data can be reliably collected.

The challenges of modern data analysis breach this standard setup. A first difference, one that cannot be overstated, is that data are being collected everywhere, on a more massive scale than ever before, by cameras, sensors, and people. We give just a few examples: There are an estimated minimum 10,000 surveillance cameras in the city of Chicago and an estimated 500,000 in London . Netflix collects ratings from 25 million users on tens of thousands of movies . On its peak day of the holiday season in 2008, Amazon.com collected data on 72 items purchased every second . The Large Synoptic Survey Telescope, which will be deployed in Chile and will photograph the whole sky visible to it every three nights, will produce 20 terabytes of data every night .

A second and equally important difference is that, in all these examples mentioned, the data collected may be unreliable or an indirect indicator of what one really wants to know. The data are collected from many possibly distributed sensors or even from people whose responses may be inconsistent, and the data may be missing or corrupted.

In order to address these issues, algorithms for data analysis must be computationally fast as well as robust to corruption and missing data. In this paper we present the Grassmannian Robust Adaptive Subspace Tracking Algorithm, or GRASTA, an online algorithm for robust subspace tracking that handles these three challenges at once. We seek a low-rank model for data that may be corrupted by outliers and have missing data values.

GRASTA uses the natural $l^{1}$-norm cost function for data corrupted by sparse outliers, and performs incremental gradient descent on the Grassmannian, the manifold of all $d$-dimensional subspaces for fixed $d$. For each subspace update, we use the gradient of the augmented Lagrangian function associated to this cost. GRASTA operates only one data vector at a time, making it faster than other state-of-the-art algorithms and amenable to streaming and real-time applications.

The contributions of our work are threefold:

We propose an efficient online robust subspace tracking algorithm – GRASTA, or Grassmannian Robust Adaptive Subspace Tracking Algorithm – which combines the augmented Lagrangian function with the classic stochastic gradient framework and the structure of the Grassmannian , and solves via the augmented Lagrangian alternating direction method . As we discuss in detail in Section 2 and 3, GRASTA alternates between estimating a low-dimensional subspace ${\mathcal{S}}$ and a triple $(s,w,y)$ which represent the sparse corruptions in the signal, the weights for the fit of the signal to the subspace, and the dual vector in the optimization problem. For estimating the subspace ${\mathcal{S}}$, GRASTA uses gradient descent on the Grassmannian with $(s,w,y)$ fixed; for estimating the triple $(s,w,y)$, GRASTA uses ADMM .

When data vectors arise from an underlying subspace which is inherently low-dimensional, and are corrupted with noise and outliers, GRASTA is able estimate and track the subspace successfully, even when the vectors are highly incomplete.

1.2 Fast Robust Low-rank Matrix Completion

We show that GRASTA can successfully recover a low-rank matrix from partial information, even if the partially observed entries are corrupted by gross outliers. GRASTA’s incremental update results in a significant speed-up over other state-of-the-art robust matrix completion algorithms or RPCA (robust principal components analysis) algorithms.

1.3 Realtime Separation of Background and Moving Objects in Video Surveillance

Finally, we show that the online nature of GRASTA makes it suitable for realtime high-dimensional sparse signal separation from a background signal, such as the task of separating background and moving objects in video surveillance. Compared to other RPCA methods, GRASTA can handle video frames at very high rates– up to 57 frames per second in our examples– even when implemented in MATLAB on a personal laptop, which is a significant practical advantage over other state-of-the-art techniques.

This paper is organized as follows. We motivate robust online subspace tracking and give background on subspace tracking and matrix completion in Sections 1.2 and 1.3. The familiar reader can go directly to Section 2, where we formulate the robust subspace tracking problem and introduce the novel subspace error function in Section 2. In Section 3, we present the Grassmannian Robust Adaptive Subspace Tracking Algorithm (GRASTA) in detail and discuss critical parts of the implementation; we point out the limitations and merits as compared with other RPCA algorithms. In Section 4, we compare GRASTA with GROUSE and RPCA algorithms via extensive numerical experiments and several real-world video surveillance experiments. Section 5 concludes our work and gives some discussion on future directions.

2 Motivations

GRASTA is built on GROUSE , an efficient online subspace tracking algorithm. GROUSE uses an $l^{2}$-norm cost function, which is problematic when facing data corruption or noise distributed other than Gaussian.

As an example, we consider using subspaces to detect anomalies in computer networks . A non-robust subspace estimation algorithm like GROUSE would need a special anomaly detection component in order to differentiate anomalies and outliers from the underlying subspace of the traffic data. Often these types of anomaly detection components rely on a lot of parameter tuning and heuristic rules for detection. This motivates a more principled approach that is robust by design: GRASTA.

2.2 Robust Principal Component Analysis

Principal Components Analysis is a critical tool for data analysis in many fields. Given a parameter $d$ for the number of components desired, PCA seeks to find the best-fit (in an $l^{2}$ norm sense) $d$-dimensional subspace to data; in other words, it finds the best $d$ vectors, the principal components, such that the data can be approximated by a linear combination of those $d$ vectors.

The residuals of an $l^{2}$-norm error function will be Gaussian distributed. Therefore, even with one outlier data point, the principal components can be arbitrarily far from those without the outlier data point . Modern data applications– such as those in sensor networks, collaborative filtering, video surveillance or the network monitoring example just given– will all experience data failures that result in outliers. Sometimes the outliers are even the signal of interest, as in the case of network anomaly detection or identifying moving objects in the foreground of a surveillance camera.

A good deal of research is therefore focused on Robust PCA, including . Recent work focuses on a problem definition which seeks a low-rank and sparse matrix whose sum is the observed data. The majority of algorithms use SVD (singular value decomposition) computations to perform Robust PCA. The SVD is too slow for many real-time applications, and consequently many online SVD and subspace identification algorithms have been developed, as we discuss in Section 1.3.1. We are therefore motivated to bridge the gap between online algorithms and robust algorithms with GRASTA.

Of course we emphasize that besides the ability to do matrix separation into low-rank and sparse parts, GRASTA can also effectively handle the scenario where the low-rank subspace is dynamic.

3 Background

First we briefly describe the subspace tracking problem set-up and GROUSE algorithm before reviewing previous literature on subspace tracking.

Comon and Golub give an early survey of adaptive methods for tracking subspaces, both coming from the matrix computation literature, including Lanczos-based recursion algorithms, and gradient-based methods from the signal processing literature.

There is a vast literature on the adaptation of QR and SVD factorizations to the adaptive, online context. The work in are all along these lines. The fastest algorithm for incremental SVD is given in ; this algorithm makes modifications, one column at a time, to the thin SVD of a strictly rank-$d$ $n\times n$ matrix in $O(n^{2}d)$ time.

Initial work in signal processing for subspace tracking was aimed at estimating from data the largest eigensubspace for a signal covariance matrix. This is useful, for example, in direction-of-arrival (DOA) estimation: the well-known work in introduces ESPRIT, a parameter estimation algorithm that estimates the DOA of plane waves emanating from a target and being received by a sensor array. ESPRIT was a follow up to the MUSIC algorithm , and ESPRIT gains computational efficiency over MUSIC for a slight tradeoff in generality of sensor array design. Around the same time, Yang and Kaveh introduced an approach for subspace tracking that, like GROUSE, uses incremental gradient, thus making it more suitable for adaptive estimation of the signal subspace and covariance matrix. This work was followed by with various improvements and convergence analyses. Unlike GROUSE, these algorithms all conduct gradient descent in the ambient space as opposed to operating along the geodesics of the Grassmannian. Also unlike GROUSE, these algorithms all require fully observed vectors.

Smith thoroughly pursued conjugate gradient descent methods on the Grassmannian for solving the subspace tracking problem using the Rayleigh quotient as a cost function as opposed to the Frobenius norm of GROUSE. In the authors give a very careful definition of the problem, giving a nice survey comparing the applicability of various approaches. In is an extensive list of subspace tracking references.

We note here that none of the work in this subsection addressed issues of robustness to corrupted data or missing data.

The work of addresses the problem of robust online subspace tracking. They focus on the problem where outliers are found in a fraction of vectors (that is, some vectors have no outliers), though they do remark that this can be extended to handle the case where outliers are sparse in every vector. They have a very nice proposition relating $l^{0}$-(pseudo)norm minimization to the least trimmed squares estimator.

We note here that GRASTA differs from in that it directly focuses on the case where every vector may have outliers, it operates on the Grassmannian for greater efficiency, and it can handle missing data. A comparison to is a subject of future investigation.

3.2 Matrix Completion

The popular Netflix prize stimulated research on the matrix completion problem: Given very few entries of a low-rank matrix, can one recover (or complete) the entire matrix? When Candès and Recht proved that, under some incoherence conditions, nuclear norm minimization recovers a highly incomplete low rank matrix with high probability , an entire area was opened up for further analysis and algorithmic variations. Algorithms that have been proposed to solve matrix completion include ADMiRA , OptSpace , Singular Value Thresholding , FPCA , SET , APGL , GROUSE , and many others. Of these, GROUSE is the only online matrix completion algorithm in that it proceeds incrementally, one column at a time. This along with the fact that each update of GROUSE has low computational complexity makes GROUSE the fastest of the state-of-the-art matrix completion algorithms by nearly an order of magnitude .

Problem Set-up

At each time step $t$, we assume that $v_{t}$ is generated by the following model:

where $w_{t}$ is the $d\times 1$ weight vector, $s_{t}$ is the $n\times 1$ sparse outlier vector whose nonzero entries may be arbitrarily large, and $\zeta_{t}$ is the $n\times 1$ zero-mean Gaussian white noise vector with small variance. We observe only a small subset of entries of $v_{t}$, denoted by $\Omega_{t}\subset\{1,\dots,n\}$.

It was shown in that this cost function gives an accurate estimate of the same cost function with full data ($\Omega=\{1,\dots,n\}$), as long as $|\Omega_{t}|$ is large enough In the authors show that $|\Omega_{t}|$ must be larger than $\mu({\mathcal{S}})d\log(2d/\delta)$, where $\mu({\mathcal{S}})$ is a measure of incoherence on the subspace and $\delta$ controls the probability of the result. See the paper for details.. However, if the observed data vector is corrupted by outliers as in Equation (2.1), an $l^{2}$-based best-fit to the subspace can be influenced arbitrarily with just one large outlier; this in turn will lead to an incorrect subspace update in the GROUSE algorithm, as we demonstrate in Section 4.1.

2 Subspace Error Quantification by l1l^{1}-Norm

In order to quantify the subspace error robustly, we use the $l^{1}$-norm as follows:

With $U_{\Omega_{t}}$ known (or estimated, but fixed), this $l^{1}$ minimization problem is the classic least absolute deviations problem; Boyd has a nice survey of algorithms to solve this problem and describes in detail a fast solver based on the technique of ADMM (Alternating Direction Method of Multipliers) http://www.stanford.edu/~boyd/papers/admm/. More references can be found therein.

According to , we can rewrite the right hand of Equation (2.3) as the equivalent constrained problem by introducing a sparse outlier vector $s$:

The augmented Lagrangian of this constrained minimization problem is then

where $y$ is the dual vector. Our unknowns are $s$, $y$, $U$, and $w$. Note that since $U$ is constrained to a non-convex manifold ($U^{T}U=I$), this function is not convex (neither is Equation (2.2)). However, note that if $U$ were estimated, we could solve for the triple $(s,w,y)$ using ADMM; also if $(s,w,y)$ were estimated, we could refine our estimate of $U$. This is the alternating approach we take with GRASTA. We describe the two parts in detail in Sections 3.1 and 3.2.

3 Relation to Robust PCA and Robust Matrix Completion

The global version of the $l^{1}$ cost function in Equation (2.3) follows:

The right hand of Equation (2.6) can be rewritten as the equivalent constrained problem:

which is the same problem studied in , and the authors propose an efficient ADMM solver for this problem. Unlike the set-up of , this problem is not convex; however it offers much more computationally efficient solutions. GRASTA differs from the algorithm of in two major ways: it uses incremental gradient to minimize this cost function one column at a time for even greater efficiency, and it uses geodesics on the Grassmannian to compute the update of $U$.

Grassmannian Robust Adaptive Subspace Tracking

As we have said, GRASTA alternates between estimating the triple $(s,w,y)$ and the subspace $U$. Here we discuss those two pieces of our algorithm. Section 3.1 describes the update of $(s,w,y)$ based on an estimate $\widehat{U}_{t}$ for the subspace variable. Section 3.2 describes the update of our subspace variable to $\widehat{U}_{t+1}$ based on the estimate of $(s^{*},w^{*},y^{*})$ resulting from the first step. Finally, Section 3.4 describes our algorithm for adaptively choosing the gradient step-size.

Given the current estimated subspace $\widehat{U}_{t}$, the partial observation $v_{\Omega_{t}}$, and the observed entries’ indices $\Omega_{t}$, the optimal $(s^{*},w^{*},y^{*})$ of Equation (2.4) can be found with the following minimization of the augmented Lagrangian.

Equation (3.1) can be efficiently solved by ADMM . That is, $s$, $w$, and the dual vector $y$ are updated in an alternating fashion:

Specifically, these quantities are computed as follows. In this paper we always assume that $U_{\Omega_{t}}^{T}U_{\Omega_{t}}$ is invertible, which is guaranteed if $|\Omega_{t}|$ is large enough . We have:

where $\mathsf{S}_{\frac{1}{1+\rho}}$ is the elementwise soft thresholding operator . We discuss this ADMM solver in detail as Algorithm 2 in Section 3.5.

2 Subspace Update

GRASTA achieves online robust subspace tracking by performing incremental gradient descent on the Grassmannian step by step. That is, we first compute a gradient of the loss function, and then follow this gradient along a short geodesic curve on the Grassmannian. Figure 1 illustrates the basic idea of gradient descent along a geodesic.

It seems that it would be natural to use Equation (2.3) as the robust loss function. However, there is a critical limitation of this approach: when regarding $U$ as the variable, this loss function is not differentiable everywhere.

Here we propose to use the augmented Lagrangian as the subspace loss function once we have estimated $(s^{*},w^{*},y^{*})$ from the previous $\widehat{U}_{\Omega_{t}}$ and $v_{\Omega_{t}}$ by Equation (3.2). The new loss function is stated as Equation (3.6):

This new subspace loss function is differentiable. Furthermore, when the data vector is not corrupted by outliers, Equation (3.6) reduces to the $l^{2}$-norm loss function of GROUSE .

2.2 Grassmannian Geodesic Gradient Step

Then the gradient $\triangledown{\mathcal{L}}$ is $\triangledown{\mathcal{L}}=(I-UU^{T})\frac{d\mathcal{L}}{dU}$ . Here we introduce three variables $\Gamma$, $\Gamma_{1}$, and $\Gamma_{2}$ to simplify the gradient expression:

Thus the gradient $\triangledown{\mathcal{L}}$ can be further simplified to:

From Equation (3.11), it is easy to verify that $\triangledown{\mathcal{L}}$ is rank one since $\Gamma$ is a $n\times 1$ vector and $w^{*}$ is the optimal $d\times 1$ weight vector. Then it is trivial to compute the singular value decomposition of $\triangledown{\mathcal{L}}$, which will be used for the following gradient descent step along the geodesic according to Equation (2.65) in . The sole non-zero singular value is $\sigma=\|\Gamma\|\|w^{*}\|$, and the corresponding left and right singular vectors are $\frac{\Gamma}{\|\Gamma\|}$ and $\frac{w^{*}}{\|w^{*}\|}$ respectively. Then we can write the SVD of the gradient explicitly by adding the orthonormal set $x_{2},\ldots,x_{d}$ orthogonal to $\Gamma$ as left singular vectors and the orthonormal set $y_{2},\ldots,y_{d}$ orthogonal to $w^{*}$ as right singular vectors as follows:

Finally, following Equation (2.65) in , a gradient step of length $\eta$ in the direction $-\triangledown{\mathcal{L}}$ is given by

3 Remarks

Here we point out that at each subspace update step, our approach does not remove outliers explicitly. In fact, we use the gradient of the augmented Lagrangian $\mathcal{L}(U)$ Equation (3.6) which exploits the dual vector $y^{*}$ to leverage the outlier effect. That is the key to success. Even when the ADMM solver 3.2 can not identify the outliers due to our current estimated subspace being far away from the true subspace, with the help of the dual vector $y^{*}$ the gradient of the augmented Lagrangian gives us the right direction at each step which leads us to the right subspace.

We also must point out that since we estimate $(s^{*},w^{*},y^{*})$ at each step using the ADMM solver, we can not recover the exact subspace with sufficient accuracy if we do not allocate enough iterations for the ADMM solver . Fortunately, as it also emphasized in , only a few tens of iterations per subspace update step are sufficient to achieve a modest accuracy, which is often acceptable for practical use. Extensive experiments in Section 4 show that our algorithm is fast and always produces acceptable results, even when the vectors are noisy and heavily corrupted by outliers.

4 Adaptive Step-size

The question of how large a gradient step to take along the geodesic is an important issue, and it depends on a fundamental tradeoff between tracking rate and steady-state error. Rather than the constant step-size rule proposed for subspace tracking in GROUSE, here we propose to use the adaptive step-size rule to achieve both precise convergence for a stationary subspace and fast adaptation to a changing subspace.

We use the following formula to update the step-size $\eta_{t}$:

where $C$ is the predefined constant step-size scale. If we use $\mu_{t}=t$ to update $\eta_{t}$, it is obvious that the step-size satisfies the following properties:

This is the classic diminishing step-size rule in stochastic gradient descent literature, and has been proven to guarantee convergence to a stationary point .

However, our goal is not only to identify the stationary subspace precisely. We have the more ambitious goal of keeping track of the subspace when the subspace is slowly changing. Obviously, with a changing subspace, if we use a diminishing step-size rule, when $\eta_{t}$ is shrinking to $0$ our steps will be too small to track the dynamic subspace. To continually adapt to the changing subspace, GROUSE proposes a constant step-size which needs careful selection to balance the tradeoff between tracking rate and steady-state error.

Here we propose to use an adaptive step-size rule to produce a proper step-size $\eta_{t}$ that empirically achieves both precise convergence for a stationary subspace and fast adaptation to a changing subspace. The basic idea is inspired by Plakhov and Klein : if two consecutive gradients $\triangledown\mathcal{L}_{t-1}$ and $\triangledown\mathcal{L}_{t}$ are in the same direction, i.e. $\langle\triangledown\mathcal{L}_{t-1},\triangledown\mathcal{L}_{t}\rangle>0$, it intuitively means that the current estimated $\widehat{U}_{t}$ is relatively far away from the true subspace $\mathcal{S}_{t}$. If this is the case, heuristically we should take a slightly larger step along $-\triangledown\mathcal{L}_{t}$ than the previous step-size $\eta_{t-1}$. Otherwise, if $\triangledown\mathcal{L}_{t-1}$ and $\triangledown\mathcal{L}_{t}$ are not in the same direction, i.e. $\langle\triangledown\mathcal{L}_{t-1},\triangledown\mathcal{L}_{t}\rangle<0$, again intuitively this means that the current estimated $\widehat{U}_{t}$ is relatively close the true subspace $\mathcal{S}_{t}$, and again heuristically we should take a slightly smaller step along $-\triangledown\mathcal{L}_{t}$ than the previous step-size $\eta_{t-1}$. Besides the sign of the two consecutive gradients giving us intuition for the step-size adaptation, the inner product $\langle\triangledown\mathcal{L}_{t-1},\triangledown\mathcal{L}_{t}\rangle$ also gives us the proper adapted magnitude for our step-size .

We still use Equation (3.14) to produce $\eta_{t}$ at each time, but update $\mu_{t}$ according to the inner product of two consecutive gradients $\langle\triangledown\mathcal{L}_{t-1},\triangledown\mathcal{L}_{t}\rangle$ as follows:

where the $sigmoid$ function is defined as:

with $sigmoid(0)=0$, $f_{MAX}>0$, $f_{MIN}<0$, and $\omega>0$. $f_{MAX}$ and $f_{MIN}$ are chosen to control how much the step-size grows or shrinks; and $\omega$ controls the shape of the $sigmoid$ function. In this paper we always set $f_{MAX}=1$, $f_{MIN}=-1$, and $\omega=0.1$.

When the estimated subspace $\widehat{U}_{t}$ is very close to the true subspace ${\mathcal{S}}_{t}$, the adaptive step-size $\eta_{t}\rightarrow 0$, or equivalently $\mu_{t}\rightarrow+\infty$ from Equation (3.14). Now we consider the following scenario: suppose we have identified the subspace precisely– and therefore $\mu_{t}>N$ for some large number $N$ then suddenly the subspace changes dramatically. How quickly will this step-size rule adapt to the new subspace? In practical applications, taking too much time to adapt to the new subspace is undesirable. Specifically, only shrinking $\mu_{t}$ at most $|f_{MIN}|$ is too conservative in this scenario. It is easy to verify that, since at each update step $\mu_{t}$ shrinks at most $|f_{MIN}|$, the increase of $\eta_{t}$ is limited and therefore this approach wouldn’t take very large steps even though the subspace has changed. When the subspace changes drastically, we should shrink $\mu_{t}$ more to accelerate the adaptation process.

For GRASTA, we take this approach and call it a ”Multi-Level” adaptive step-size rule. Though we do not provide the convergence proof here, empirically this multi-level adaptive approach demonstrates much faster convergence performance than the single-level strategy discussed above. We leave further detailed comparison to future investigation.

Our multi-level adaptation is as follows. We only let $\mu_{t}$ change in $\left(\mu_{MIN},\mu_{MAX}\right)$, where $\mu_{MIN}$ and $\mu_{MAX}$ are prescribed constants. For the experiments in this paper we always set $\mu_{MIN}=1$ and $\mu_{MAX}=15$. Then in this case Equation (3.15) is adapted to Equation (3.16):

We introduce a level variable $l_{t}$ that will get smaller when our subspace estimate is far from the data. Then the step-size $\eta_{t}$ is as follows:

Once $\mu_{t}$ calculated by Equation (3.16) is larger than $\mu_{MAX}$, we increase the level variable $l_{t}$ by $1$ and set $\mu_{t}=\mu_{0}$, where $\mu_{0}\in\left(\mu_{MIN},\mu_{MAX}\right)$ and $\mu_{0}$ is selected close to $\mu_{MIN}$ (in our experiments we let $\mu_{0}=3$). If $\mu_{t}\leq\mu_{MIN}$, we decrease $l_{t}$ by $1$ and also set $\mu_{t}=\mu_{0}$. Therefore, when our subspace estimate is off, we are increasing $\eta_{t}$ exponentially instead of linearly. On one hand this new multi-level adaptive rule follows the basic adaptive step-size rule discussed above; on the other hand exploiting this multi-level property, this new approach adapts more quickly to a changing subspace. Once we have identified the subspace changing and $\mu_{t}\leq\mu_{MIN}$, if the subspace really changes dramatically, $l_{t}$ will keep decreasing until $\mu_{t}$ is again within the range $\left(\mu_{MIN},\mu_{MAX}\right)$.

Combining these ideas together, we state our novel adaptive step-size rule as Algorithm 3.

5 Algorithms

The discussion of Sections 3.1 to 3.4 can be summarized into our algorithm as follows. For each time step $t$, when we observe an incomplete and corrupted data vector $v_{\Omega_{t}}$, our algorithm will first estimate the optimal value $(s^{*},w^{*},y^{*})$ from our current estimated subspace $U_{t}$ via the $l^{1}$ minimization ADMM solver 3.2; then compute the gradient of the augmented Lagrangian loss function $\mathcal{L}$ by Equation (3.11); then estimate a proper step-size $\eta_{t}$ from the two consecutive gradients $\triangledown\mathcal{L}_{t-1}$ and $\triangledown\mathcal{L}_{t}$ by Equation (3.15) and 3.17 ; and finally do the rank one subspace update via Equation (3.13).

We state our main algorithm GRASTA (Grassmannian Robust Adaptive Subspace Tracking Algorithm) in Algorithm 1. GRASTA consists of two important sub-procedures: the ADMM solver of the least absolute derivations problem, and the computation of the adaptive step-size. We state the two sub-procedures as Algorithm 2 and Algorithm 3 separately.

Unlike GROUSE, which has a closed form solution for computing the gradient, GRASTA estimates $(s_{t}^{*},w_{t}^{*},y_{t}^{*})$ by the ADMM iterated Algorithm 2. Certainly we would have a potential performance bottleneck if Algorithm 2 takes too much time at each subspace update step. However, we see empirically that only a few tens of iterations in Algorithm 2 at each step allows GRASTA to track the subspace to an acceptable accuracy. In our video experiments with Algorithm 2, we always set the maximum iteration $K$ around $20$ to balance the trade-off between the subspace tracking accuracy and computational performance. We make a slight modification to the original ADMM sovler presented in : in addition to returning $w^{*}$ we also return the sparse vector $s^{*}$ and the dual vector $y^{*}$ for the further computation of the gradient $\triangledown{\mathcal{L}}$. It is easy to verify that in the worst case the ADMM solver needs at most $O(|\Omega|d^{3}+Kd|\Omega|)$ flops.

In order to produce the proper step-size $\eta_{t}$ from Algorithm 3, we need to maintain the gradient $\triangledown\mathcal{L}_{t-1}$ from the previous time step throughout the subspace tracking process. Keeping $\triangledown\mathcal{L}_{t-1}$ only requires additional $O(n+d)$ memory usage. The main computation of of Algorithm 3 is the inner product $\langle\triangledown\mathcal{L}_{t-1},\triangledown\mathcal{L}_{t}\rangle$, which is the trace of the product $\triangledown\mathcal{L}_{t-1}$ and $\triangledown\mathcal{L}_{t}$, two $n\times d$ matrices, and will cost $O(nd^{2})$ flops.

6 Computational Cost and Memory Usage

Each subspace update step in GRASTA needs only simple linear algebraic computations. The total computational cost of each step of Algorithm 1 is $O(|\Omega|d^{3}+Kd|\Omega|+nd^{2})$, where again $|\Omega|$ is the number of samples per vector used, $d$ is the dimension of the subspace, $n$ is the ambient dimension, and $K$ is the number of ADMM iterations.

Specifically, estimating $(s_{t}^{*},w_{t}^{*},y_{t}^{*})$ from Algorithm 2 costs at most $O(|\Omega|d^{3}+Kd|\Omega|)$ flops; computing the gradient $\triangledown\mathcal{L}$ needs simple matrix-vector multiplication which costs $O(|\Omega|d+nd)$ flops; producing the adaptive step-size costs $O(nd^{2})$ flops; and the final update step also costs $O(nd^{2})$ flops.

Throughout the tracking process, GRASTA only needs $O(nd)$ memory elements to maintain the estimated low-rank orthonormal basis $\widehat{U}_{t}$, $O(n)$ elements for $s^{*}$ and $y^{*}$, $O(d)$ elements for $w^{*}$, and $O(n+d)$ for the previous step gradient $\triangledown\mathcal{L}_{t-1}$ in memory.

This analysis decidedly shows that GRASTA is both computation and memory efficient.

Numerical Experiments

In the following experiments, we explore GRASTA’s performance in various scenarios: subspace tracking, robust matrix completion, and the video surveillance application. We use relative error to quantify the performance of GRASTA. If the recovered data is a vector, the relative error is defined as follows:

If the recovered data is a matrix, the relative error is defined as follows:

We also use ”Noise Relative Power” to quantify the additional Gaussian white noise perturbation, which is defined as follows:

Here $v$ is the true data vector and $\zeta$ is the additional Gaussian noise as in Equation (2.1).

In all the following experiments, we use Matlab R2010b on a Macbook Pro laptop with 2.3GHz Intel Core i5 CPU and 4 GB RAM. To improve the performance, we implement Algorithm 2 in C++ and make it as a MEX-file to be integrated into GRASTA Matlab scripts.

Our first goal is to compare GRASTA with the non-robust algorithm GROUSE to show the need for a robust subspace estimation and tracking algorithm.

In many of the following experiments, we use this generative model to generate a series of data vectors:

$U_{true}$ is an $n\times d$ matrix whose $d$ columns are realizations of an i.i.d. $\mathcal{N}(0,I_{n})$ random variable that are then orthornomalized. The weight vector $w_{t}$ is a $d\times 1$ vector whose entries are realizations of i.i.d. $\mathcal{N}(0,1)$ random variables, that is Gaussian distributed with mean zero and variance 1. The sparse vector $s_{t}$ is an $n\times 1$ vector whose nonzero entries are Gaussian noise with the maximum of the data vector $U_{true}w$ as the variance; the locations of the nonzero entries are chosen uniformly at random without replacement. The noise $\zeta_{t}$ is an $n\times 1$ vector whose entries are i.i.d $\mathcal{N}(0,\omega^{2})$. This parameter $\omega^{2}$ governs the SNR with respect to the low-rank part of our data. For the entire comparison against GROUSE, we used a maximum of $K=60$ iterations of the ADMM algorithm per subspace iteration.

Figure 2 illustrates the failure of GROUSE, and success of GRASTA, when these sparse outliers are added only at periodic time intervals. We can see that GROUSE is significantly thrown off, despite the outliers occurring in an isolated vector. This illustrates clearly our motivation for adding robustness to the subspace tracking algorithm.

1.2 Robust Matrix Completion

We aim to complete $500\times 500$ dimensional matrices of rank $5$. The matrices are corrupted by different fractions of outliers, depending on the experimental setting, and sampled uniformly without replacement with density $0.30$. We generate the low-rank matrix by first generating two $500\times 5$ factors $Y_{L}$ and $Y_{R}$ with i.i.d. Gaussian entries and then adding normally distributed noise with variance $\omega^{2}$. The location of sparse outliers is distributed uniformly, and the outlier values are normally distributed with variance equal to the maximum of the matrix.

For each setting of the fraction of outliers, we randomly generate $5$ matrices, each of which is solved via GROUSE and GRASTA separately. Both GROUSE and GRASTA cycle through the matrix columns $10$ times. Table 1 shows the averaged results of a comparison between GROUSE and GRASTA. As expected, GRASTA vastly outperforms GROUSE across the board even with the smallest number of outliers.

2 Stationary Subspace Identification

Now we wish to examine GRASTA’s performance on the stationary subspace identification problem under various conditions. In most experiments (and unless otherwise noted) the ambient dimension is $n=500$ and the inherent subspace dimension is $d=5$. We again generate the vectors using Equation (4.4) above and the descriptive text that follows Equation (4.4). We vary the fraction of entries that are corrupted, and we vary the fraction of entries that are observed.

We start with Figure 3, which shows subspace estimation performance under a varying fraction of added outliers. We can see in this problem instance that with 10% corrupted entries, the relative error reaches the relative noise floor after a number of iterations that is a small multiple of the ambient dimension. For more corruption, more vectors (gradient iterations) are needed, but even with 50% outliers and more, the relative error trends toward the relative noise power.

In Figure 4, we consider GRASTA’s error performance for varying sub-sampling rates. Here the fraction of corrupted values is fixed at 10%. We can see that again, even with a 30% sampling rate, the relative error quickly reaches the relative noise power.

Now we wish to take a closer look at the case when we have both dense outlier corruption and subsampling of the signal. This is an important scenario for applications where the “outlier corruption” is a signal of interest obscuring a low-rank background signal, and we wish to subsample in order to improve computational complexity. For example this would apply to anomaly detection problems or to the problem of separation of background and moving objects in video as we show in Section 4.5.

Figure 5 illustrates that even when the vector is highly corrupted with 50% added outliers, GRASTA can identify the underlying low-rank subspace even with only 50% of the entries. We vary the dimension (or rank) of the underlying subspace, and because of this there is not one relative noise power benchmark to compare against; however we see that the trend is similar to those in previous figures.

3 Dynamic Subspace Tracking

The fact that GRASTA operates one vector at a time allows it to track an evolving subspace. In this section we show GRASTA’s performance under two models of evolving subspaces. In these experiments, we use the same set-up as before: $n=500$, $d=5$, and $v_{t}$ is generated by Equation (4.4), except that $U_{true}=U[t]$, i.e. the subspace we wish to estimate varies with time $t$:

We use the following ordinary differential equation to simulate a rotating subspace:

where $B$ is a skew-symmetric matrix. Consequently, the subspace $U[t]$ is updated via

where $\delta$ controls the amount of rotation of with each time step $t$. As we see in Figure 6, for the rotation parameter $\delta$ fixed at $10^{-5}$, GRASTA successfully latches on and tracks the rotating subspace.

3.2 Sudden Subspace Change Tracking

For this experiment, we wanted to see the behavior of GRASTA when the subspace experienced a sudden dramatic change. At intervals of 5000 vectors, we randomly changed the true subspace to a new random subspace. The results are in Figure 7. Again from these simulations we see that GRASTA successfully identifies the subspace change and tracks the subspace again.

4 Comparison with Robust PCA

Here we compare GRASTA with RPCA on the recovery of corrupted low-rank matrices. For RPCA we use , or the IALM (Inexact Augmented Lagrange Multiplier) method The code we used is available here: http://perception.csl.uiuc.edu/matrix-rank/sample_code.html. We downloaded it in April 2011..

The corrupted matrices can be written as $M=L+S+N$, where $L$ are the low-rank matrices we want to recover, $S$ are the sparse outlier matrices, and $N$ are the Gaussian noise matrices with small variance relative to the sparse outliers. We use $d=5$ matrices of size $2000\times 2000$ to do the comparison. The low-rank matrices $L$ are generated by the same method of the previous robust matrix completion experiments: as the product of two $2000\times 5$ factors $Y_{L}$ and $Y_{R}$ with i.i.d. Gaussian entries. The sparse outlier matrices $S$ are generated by selecting a fraction of entries uniformly at random without replacement, whose values are set according to Gaussian distribution with the maximum of $L$ as the variance We note here that in , the authors use a uniform distribution for the outliers, as opposed to Gaussian. The authors in use $\pm 1$ Bernoulli variables. Gaussian is the most challenging case, because more outliers will be near zero and confuse the estimation.. We vary the fraction of corruptions from $10\%$ (sparse outliers) to $50\%$ (dense outliers), and we also vary the variance Gaussian noise matrices $N$ from a moderate perturbation of $\omega^{2}=10^{-4}$ to a larger perturbation of $\omega^{2}=10^{-3}$. For GRASTA we cycled through the matrix columns twice and used a maximum of $K=60$ iterations of the ADMM algorithm; we used a maximum of $20$ iterations of IALM.

Table 2 shows the results of the comparison. We ran RPCA with full data, GRASTA with full data, and then GRASTA with various levels of subsampling. When there are very few outliers and little noise, RPCA achieves a reasonable error rate at computational speeds similar to GRASTA. However with an increase in noise or fraction of outliers, GRASTA achieves good error performance in much less time. As a particular example, when $\omega^{2}=10^{-3}$ and the fraction of outliers is 30%, in 69 seconds GRASTA with full data achieves better error performance than RPCA in 363 seconds, and GRASTA with 30% subsampling achieves better error performance in only 23 seconds.

5 Realtime Video Background Tracking and Foreground Detection

In this subsection we discuss the application of GRASTA to the prominent problem of realtime separation of foreground objects from the background in video surveillance. Imagine we had a video with only the background: When the columns of a single frame of this video are stacked into a single column, several frames together will lie in a low-dimensional subspace. In fact if the background is completely static, the subspace would be one-dimensional. That subspace can be estimated in order to identify and separate the foreground objects; if the background is dynamic, subspace tracking is necessary. GRASTA is uniquely suited for this burgeoning application.

Here we consider three scenarios in the video tasks, with a spectrum of challenges for subspace tracking. In the first we have a video with a static background and objects moving in the foreground. In the second, we have a video with a still background but with changing lighting. In the third, we simulate a panning camera to examine GRASTA’s performance with a dynamic background. The results are summarized in Table 3.

If the video background is known to be static or near static, we can use GRASTA to track the background and separate the moving foreground objects in real-time. Since the background is static, we use GRASTA first to identify the background, and then we use only Algorithm 2 to separate the foreground from the background. More precisely we do the following:

Randomly select a few frames of the video to train the static low-rank subspace $U$. In our experiments, we select frames randomly from the entire video; however for real-time processing these frames may be chosen from initial piece of the video, as long as we can be confident that every pixel of the background is visible in one of the selected frames. The low-rank subspace $U$ is then identified from these frames using partial information. We use $30\%$ of the pixels, select $50$ frames for training, and set RANK = 5 in all the following experiments.

Once the video background $BG$ has been identified as a subspace $U$, separating the foreground objects $FG$ from each frame can be simply done using Equation (4.7), where the weight vector $w_{t}$ can be solved for via Algorithm 2, again from a small subsample of each frame’s pixels.

Table 3 shows the real-time We comment here that to call something “real-time” processing of course will depend on one’s application requirements and hardware (camera frame capture rate, in the example of video processing). For example, standard 35mm film video uses 24 unique frames per second. The maximum frame rate for most CCTVs is 30 frames per second. video separation results. From the first experiment, we use the “Hall” dataset from which consists of $3584$ frames each with resolution $144\times 176$. We let GRASTA cycle $5$ times over the $50$ training frames just from $30\%$ random entries of each frame to get the stationary subspace $U$. Training the subspace costs $6.9$ seconds. Then we perform background and foreground separation for all frames in a streaming fashion, and when dealing with each frame we only randomly observe $5\%$ entries. The separation task is performed by Equation (4.7), and the separating time is $62.5$ seconds, which means we achieve $57.3$ FPS (frames per second) real-time performance. Figure 8 shows the separation quality at $t=1,230,1400$. In order to show GRASTA can handle higher resolution video effectively, we use the “Shopping Mall” video with resolution $320\times 256$ as the second experiment. We also do the subspace training stage with the same parameter settings as “Hall”. We do the background and foreground separation only from $1\%$ entries of each frame. For “Shopping Mall” the separating time is $39.1$ seconds for total $1286$ frames. Thus we achieve $32.9$ FPS real-time performance. Figure 9 shows the separation quality at $t=1,600,1200$. In all of these video experiments we used a maximum of $K=20$ iterations of the ADMM algorithm per subspace update. The details of each tracking set-up are described in Table 4.

5.2 Dynamic Background: Changing Lighting

Here we want to consider a problem where the lighting in the video is changing throughout. We use the “Lobby” dataset from , which has $1546$ frames, each $144\times 176$ pixels. In order to adjust to the lighting changes, GRASTA tracks the subspace throughout the video; that is, unlike the last two experiments, we run the full GRASTA Algorithm 1 for every frame. We use 30% of the pixels of every frame to do this update and 100% of the pixels to do the separation. Again, see the numerical results in Table 3. The results are illustrated in Figure 10.

5.3 Dynamic Background: Virtual Pan

In the last experiment, we demonstrate that GRASTA can effectively track the right subspace in video with a dynamic background. We consider panning a ”virtual camera” from left to right and right to left through the video to simulate a dynamic background. Periodically, the virtual camera pans 20 pixels. The idea of the virtual camera is illustrated cleanly with Figure 11.

We choose “Hall” as the original dataset. The original resolution is $144\times 176$, and we set the scope of the virtual camera to have the same height but half the width, so the resolution of the virtual camera is $144\times 88$. We set the subspace $RANK=5$. Figure 12 shows how GRASTA can quickly adapt to the changed background in just $25$ frames when the virtual camera pans $20$ pixels to the right at $t=101$. We also let GRASTA track all the $3584$ frames and do the separation task for all frames. When we use 100% of the pixels for the tracking and separation, the total computation time is $191.3$ seconds, or $18.7$ FPS, and adjusting to a new camera position after the camera pans takes $25$ frames as can be seen in Figure 12. When we use 50% of the pixels for tracking and 100% of the pixels for separation, the total computation time is $144.8$ seconds or $24.8$ FPS, and the adjustment to the new camera position takes around 50 frames.

Discussion and Future Work

In this paper we have presented a robust online subspace tracking algorithm, GRASTA. The algorithm estimates a low-rank model from noisy, corrupted, and incomplete data, even when the best low-rank model may be changing over time.

Though this work presents some very successful algorithms, many questions remain. First and foremost, because the cost function in Equation (2.3) has the subspace variable $U$ which is constrained to a non-convex manifold, the resulting optimization is non-convex. A proof of convergence to the global minimum of this algorithm is of great interest.

GRASTA uses alternating minimization, alternating first to estimate $(s,w,y)$ and then fixing this triple of variables to estimate $U$. Observe that if $(s,w,y)$ are correct estimates, we could then estimate $U$ without the robust cost function. This would be quite useful in situations when speed is of utmost importance, as the GROUSE subspace update is faster than the GRASTA subspace update. Of course, knowing when $(s,w,y)$ are accurate is a very tricky business. Exploring this tradeoff is part of our future work.

We have shown that one of the very promising applications of GRASTA is that of separating background and foreground in video surveillance. We are very interested to apply GRASTA to more videos with dynamic backgrounds: for example, natural background scenery which may blow in the wind. In doing this we will study the resulting trade-off between the kinds of movement that would be captured as part of the background and the movement that would be identified as foreground.

Acknowledgments

The authors would like to thank IPAM, the Institute for Pure and Applied Mathematics, and the Internet Multi-Resolution Analysis program, which brought them together to work on this problem. We also thank Rob Nowak and Ben Recht for their thoughtful suggestions.

References