Online Identification and Tracking of Subspaces from Highly Incomplete Information

Introduction

The evolution of high-dimensional dynamical systems can often be well summarized or approximated in low-dimensional subspaces. Certain patterns of computer network traffic including origin-destination flows can be well represented by a subspace model . Environmental monitoring of soil and crop conditions , water contamination , and seismological activity have all been demonstrated to be efficiently summarized by very low-dimensional subspace representations.

The subspace methods employed in the aforementioned references are based upon first collecting full-dimensional data from the systems and then using approximation techniques to identify an accurate low-dimensional representation. However, it is often difficult or even infeasible to acquire and process full-dimensional measurements. For example, collecting network traffic measurements at a very large number of points and fine time-scales is impractical. An alternative is to randomly subsample the full-dimensional data. If the complete full-dimensional data is well-approximated by a lower dimensional subspace, and hence is in some sense redundant, then it is conceivable that the subsampled data may provide sufficient information for the recovery of that subspace. This is the central intuition of our proposed on-line algorithm for identifying and tracking low-dimensional subspaces from highly incomplete (i.e., undersampled) data. Such an algorithm could enable rapid detection of traffic spikes or intrusions in computer networks or could provide large efficiency gains in managing energy consumption in a large office building .

In this paper, we present GROUSE (Grassmannian Rank-One Update Subspace Estimation), a subspace identification and tracking algorithm that builds high quality estimates from very sparsely sampled vectors. GROUSE implements an incremental gradient procedure with computational complexity linear in dimensions of the problem, and is scalable to very high-dimensional applications. An additional feature of GROUSE is that it can be immediately adapted to an ‘online’ version of the matrix completion problem, where one aims to recover a low-rank matrix from a small, streaming random subsets of its entries. GROUSE is not only remarkably efficient for online matrix completion, but additionally enables incremental updates as columns are added or entries are incremented over time. These features are particularly attractive for maintaining databases of user preferences and collaborative filtering.

Problem Set-up

We will use this definition in Section 3 to derive our algorithm. If the matrix $U_{\Omega_{t}}^{T}U_{\Omega_{t}}$ has full rank, then we must have that $w=(U_{\Omega_{t}}^{T}U_{\Omega_{t}})^{-1}U_{\Omega_{t}}^{T}v_{\Omega_{t}}$ achieves the minimum in (1). Thus,

In the special case where the subspace is time-invariant, that is $S[t]=S_{0}$ for some fixed subspace $S_{0}$, then it is natural to consider the average cost function

The average cost function will allow us to estimate the steady-state behavior of our algorithm. Indeed, in the static case, our algorithm will be guaranteed to converge to a stationary point of $\bar{F}(S)$.

In the scenario where the subspace does not evolve over time and we only observe vectors on a finite time horizon, then the cost function (2) is identical to the matrix completion optimization problem studied in . To see the equivalence, let $\Omega=\{(k,t)~:~k\in\Omega_{t}\,1\leq t\leq T\}$, and let $V=[v_{1},\ldots,v_{T}]$. Then

That is, the global optimization problem can be written as $\min_{U,A}\sum_{(i,j)\in\Omega}(UA-V)_{ij}^{2}$, which is precisely the starting point for the algorithms and analyses in . The authors in use a gradient descent algorithm to jointly minimize both $U$ and $A$ while minimizes this cost function by first solving for $A$ and then taking a gradient step with respect to $U$. In the present work, we consider optimizing this cost function one column at a time. We show that by using our online algorithm, where each measurement $v_{t}$ corresponds to a random column of the matrix $V$, we achieve state-of-the-art performance on matrix completion problems.

Stochastic Gradient Descent on the Grassmannian

We follow the program developed in . To compute the gradient of $F$ on the Grassmannian manifold, we first need to compute the partial derivatives of $F$ with respect to the components of $U$. For a generic subspace, the matrix $U_{\Omega_{t}}^{T}U_{\Omega_{t}}$ has full rank provided that $|\Omega_{t}|>d$, and hence the cost function (1) is differentiable almost everywhere. Let $\Delta_{\Omega_{t}}$ be the $n\times n$ diagonal matrix which has $1$ in the $j^{th}$ diagonal entry if $j\in\Omega_{t}$ and has $0$ otherwise. We can rewrite

from which it follows that the derivative of $F$ with respect to the elements of $U$ is

where $r:=\Delta_{\Omega_{t}}(v_{t}-Uw)$ denotes the (zero padded) residual vector and $w$ is the least-squares solution in (1).

Using Equation (2.70) in , we can calculate the gradient on the Grassmannian from this partial derivative

The final equality follows because the residual vector $r$ is orthogonal to all of the columns of $U$. This can be verified from the definitions of $r$ and $w$.

A gradient step along the geodesic with tangent vector $-\nabla F$ is given by Equation (2.65) in , and is a function of the singular values and vectors of $\nabla F$. It is trivial to compute the singular value decomposition of $\nabla F$, as it is rank one. The sole non-zero singular value is $\sigma=2||r||||w||$ and the corresponding left and right singular vectors are $\tfrac{r}{\|r\|}$ and $\frac{w}{\|w\|}$ respectively. Let $x_{2},\ldots,x_{d}$ be an orthonormal set orthogonal to $r$ and $y_{2},\ldots,y_{d}$ be an orthonormal set orthogonal to $w$. Then

forms an SVD for the gradient. Now using (2.65) from , we find that for $\eta>0$, a step of length $\eta$ in the direction $\nabla F$ is given by

where $p:=Uw$, the predicted value of the projection of the vector $v$ onto $S$.

This geodesic update rule is remarkable for a number of reasons. First of all, it consists only of a rank-one modification of the current subspace basis $U$. Second, the term $\frac{\sin(\sigma\eta)}{\|r\|\|w\|}=\frac{\sin(\sigma\eta)}{\sigma}$ is on the order of $\eta$ when $\sigma\eta$ is small. That is, for small values of $\sigma$ and $\eta$ this expression looks like a normal step along the gradient direction $-2rw^{T}$ given by (3). From the Taylor series of the cosine, we see that the second term is approximately equal to $\sigma^{2}\eta^{2}\frac{pw^{T}}{\|p\|\|w\|}$. That is, this term serves as a second order correction to keep the iterates on the Grassmannian. Surprisingly, this simple additive term maintains the orthogonality, obviating the need for orthogonalizing the columns of $U$ after a gradient step. Below, we will also discuss how this iterate relates to more familiar iterative algorithms from linear algebra which use full information.

The GROUSE algorithm simply follows geodesics along the gradients of $F$ with a prescribed set of step-sizes $\eta$. The full computation is summarized in Algorithm 1. Our derivations have shown that computing a gradient step only requires the solution of the least squares problem (1), the computation of $p$ and $r$, and then a rank one update to the previous subspace.

Each step of GROUSE can be performed efficiently with standard linear algebra packages. Computing the weights in Step 2 of Algorithm 1 requires solving a least squares problem in $|\Omega_{t}|$ equations and $d$ unknowns. Such a system is solvable in at most $O(|\Omega_{t}|d^{2})$ flops in the worst case. Predicting the component of $v$ that lies in the current subspace requires a matrix vector multiply that can be computed in $O(nd)$ flops. Computing the residual then only requires $O(|\Omega_{t}|)$ flops, as we will always have zeros in the entries indexed by the complement of $\Omega_{t}$. Computing the norms of $r$ and $p$ can be done in $O(n)$ flops. The final subspace update consists of adding a rank one matrix to an $n\times d$ matrix and can be computed in $O(nd)$ flops. Totaling all of these computation times gives an overall complexity estimate of $O(nd+|\Omega_{t}|d^{2})$ flops per subspace update.

In the static case where the subspace $S[t]=S_{0}$ for all $t$, we can guarantee that the algorithm converges to a stationary point of (2) as long as the stepsizes $\eta_{t}$ satisfy

Selecting $\eta_{t}\propto 1/t$ will satisfy this assumption. This analysis appeals to the classical ODE method , and we are guaranteed such convergence because $\mathfrak{G}(n,d)$ is compact.

In the case that $S[t]$ is changing over time, a constant stepsize is needed to continually adapt to the changing subspace. Of course, if a non-vanishing stepsize is used then the error will not converge to zero, even in the static case. This leads to the common tradeoff between tracking rate and steady-state error in adaptive filtering problems. We explore this tradeoff in Section 4.

2 Comparison to methods that use full information

GROUSE can be understood as an adaptation of an incremental update to a QR or SVD factorization. Most batch subspace identification algorithms that rely on the eigenvalue decomposition, the singular value decomposition, or their more efficient counterparts such as the QR decomposition or the Lanczos method, can be adapted for on-line updates and tracking of the principal subspace. A comprehensive survey of these methods can be found in .

In lieu of being able to exactly able to compute component of $v_{t}$ that is orthogonal to our current subspace estimate, GROUSE computes this component only on the entries in $\Omega_{t}$. Recent work shows that for generic subspaces, the estimate for $r$ computed by Algorithm 1 in a single iteration is an excellent proxy for the amount of energy that lies in the subspace, provided that the number of measurements at each time step is greater than the true subspace dimension times a logarithmic factor. One can also verify that the GROUSE update rule corresponds to forming the matrix $[U,r]$ and then truncating the last column of the matrix

where $R_{\eta}$ denotes the $(r+1)\times(r+1)$ rotation matrix

That is, our algorithm computes a mixture of the current subspace estimate and a predicted orthogonal component. This mixture is determined both by the stepsize and the relative energy of $v_{t}$ outside the current subspace.

Numerical Experiments

We first consider the problem of identifying a fixed subspace. In all of the following experiments, the full data dimension is $n=700$, the rank of the underlying subspace is $d=10$, and the sampling density is $0.17$ unless otherwise noted. We generated a series of iid vectors $v_{t}$ according to generative model:

We implemented GROUSE (see Algorithm 1 above) with a stepsize rule of $\eta_{t}:=C/t$ for some constant $C>0$. Figure 1(a) shows the steady state error of the tracked subspace with varying values of $C$ and the noise variance. All the data points reflect the error performance at $t=14000$. We see that GROUSE converges for $C$ ranging over an order of magnitude, however with additive noise the smaller stepsizes yield smaller errors. When there is no noise, i.e., $\omega^{2}=0$, the error performance is near the level of machine precision and is flat for the whole range of converging stepsizes. In Figure 1(b) we show the number of input vectors after which the algorithm converges to an error of less than $10^{-6}$. The results are consistent with Figure 1(a), demonstrating that smaller stepsizes in the suitable range take fewer vectors until convergence. We only ran the algorithm up to time $t=14000$, so the data points for the smallest and the largest few stepsizes only reflect that the algorithm did not yet reach the desired error in the allotted time. Figures 1(c) and 1(d) repeat identical experiments with a constant stepsize policy, $\eta_{t}=C$. We again see a wide range of stepsizes for which GROUSE converges, though the region of stability is narrower in this case.

We note that the norm of the residual $\|r\|$ provides an excellent indicator for whether tracking is successful. A shown in Figure 2(a), the error to the true subspace is closely approximated by $\|r\|/\|v_{t}\|$. This confirms the theoretical analysis in which proves that this residual norm is an accurate estimator of the true subspace error when the number of samples is appropriately large.

Subspace Change Detection

As a first example of GROUSE’s ability to adapt to changes in the underlying subspace, we simulated a scenario where the underlying subspace abruptly changes at three points over the course of an experiment with 14000 observations. At each break, we selected a new subspace $S$ uniformly at random and GROUSE was implemented with a constant stepsize. As is to be expected, the algorithm is able to re-estimate the new subspace in a time depending on the magnitude of the constant stepsize.

Rotating Subspace

In this second synthetic experiment, the subspace evolves according to a random ordinary differential equation. Specifically, we sample a skew-symmetric matrix $B$ with independent, normally distributed entries and set

The resulting subspace at each iteration is thus $U[t]=\exp(\delta tB)$ where $\delta$ is a positive constant. The resulting subspace at each iteration is thus $U[t]=\exp(\delta tB)$ where $\delta$ is a positive constant. In Figure 3, we show the results of tracking the rotating subspace with $\delta=10^{-5}$. To demonstrate the effectiveness of the tracking, we display the projection of four random vectors using both the true subspace (in blue) and our subspace estimate at that time instant (in red).

Tracking Chlorine Levels

We also analyzed the performance of the GROUSE algorithm on simulated chlorine level monitoring in a pressurized water delivery system. The data were generated using EPANET http://www.epa.gov/nrmrl/wswrd/dw/epanet.html software and were previously analyzed . The input to EPANET is a network layout of water pipes and the output has many variables including the chemical levels, one of which is the chlorine level. The data we used is available from http://www.cs.cmu.edu/afs/cs/project/spirit-1/www/. This dataset has ambient dimension $n=166$, and $T=4610$ data vectors were collected, once every 5 minutes over 15 days. We tracked an $d=6$ dimensional subspace and compare this with the best 6-dimensional SVD approximation of the entire complete dataset. The results are displayed in Figure 4. The table gives the results for various constant stepsizes and various fractions of sampled data. The smallest sampling fraction we used was $20\%$, and for that the best stepsize was 3e-2; we also ran GROUSE on the full data, whose best stepsize was 5e-3. As we can see, the performance error improves for the smaller stepsize of 5e-3 as the sampling fraction increases; Also the performance error improves for the larger stepsize of 3e-2 as the sampling fraction decreases. However for all intermediate sampling fractions there are intermediate stepsizes that perform near the best reconstruction error of about 0.12. The normalized error of the full data to the best 3-dimensional SVD approximation is 0.0704. Note that we only allow for one pass over the data set and yet attain, even with very sparse sampling, comparable accuracy to a full batch SVD which has access to all of the data.

Figure 4(a) and (b) show the original and the GROUSE reconstructions of five of the chlorine sensor outputs. We plot the last 500 of the 4310 samples, each reconstructed by the estimated subspace at that time instant.

2 Matrix Completion Problems

As described in Section 2.1, matrix completion can be thought of as a subspace identification problem where we aim to identify the column space of the unknown low-rank matrix. Once the column space is identified, we can project the incomplete columns onto that subspace in order to complete the matrix. We have examined GROUSE in this context with excellent results. Our approach is to do descent on the column vectors in random order, and allow the algorithm to pass over those same incomplete columns a few times.

Our simulation set-up aimed to complete $700\times 700$ dimensional matrices of rank $10$ sampled with density $0.17$. We generated the low-rank matrix by generating two factors $Y_{L}$ and $Y_{R}$ with i.i.d. Gaussian entries, and added normally distributed noise with variance $\omega^{2}$. The robustness to step-size and time to recovery are shown in Figure 5.

In Figure 6 we show a comparison of five matrix completion algorithms and GROUSE. Namely, we compare to the performance of OPT-SPACE , FPCA , SVT , SDPLR , and NNLS . We downloaded each of these MATLAB codes from the original developer’s websites when possible. We use the same random matrix model as in Figure 5. GROUSE is faster than all other algorithms, and achieves higher quality reconstructions on many instances. We subsequently compared against NNLS, the fastest batch method, on very large problems. Both GROUSE and NNLS achieved excellent reconstruction, but GROUSE was twice as fast.

Discussion and Future Work

The inherent simplicity and empirical power of the GROUSE algorithm merit further theoretical investigations to determine exactly when it provides consistent estimates. It is possible that an adaptation of the analysis of to this online setting will yield such consistency bounds. The main difficulty lies in characterizing the trajectory of the GROUSE algorithm from a random starting subspace. Though we can guarantee that there is a basin of attraction around the global minimum, it is not yet clear how to characterize when GROUSE will end up in this basin.

We would also like to investigate how to adapt step-size in the GROUSE algorithm automatically to varying data. This is the only parameter required to run the GROUSE algorithm, and we have seen that there are substantial performance gains when this parameter is optimized. Using techniques from Least-Squares Estimation, it may be possible to automatically tune the step size based on our current error residuals.

Acknowledgments

This work was partially supported by the AFOSR grant FA9550-09-1-0140. R. Nowak would also like to thank Trinity College and the Isaac Newton Institute at the University of Cambridge for support while this work was being completed.

References