{'fraction_non_alphanumeric': 0.039405284563794545, 'fraction_numerical': 0.012120559161795998, 'mean_word_length': 4.1840268081541465, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
PETRELS: Parallel Subspace Estimation and Tracking by Recursive Least Squares from Partial Observations

Many real world data sets exhibit an embedding of low-dimensional structure in a high-dimensional manifold. Examples include images, videos and internet traffic data. It is of great significance to reduce the storage requirements and computational complexity when the data dimension is high. Therefore we consider the problem of reconstructing a data stream from a small subset of its entries, where the data is assumed to lie in a low-dimensional linear subspace, possibly corrupted by noise. We further consider tracking the change of the underlying subspace, which can be applied to applications such as video denoising, network monitoring and anomaly detection. Our problem can be viewed as a sequential low-rank matrix completion problem in which the subspace is learned in an on-line fashion. The proposed algorithm, dubbed Parallel Estimation and Tracking by REcursive Least Squares (PETRELS), first identifies the underlying low-dimensional subspace via a recursive procedure for each row of the subspace matrix in parallel with discounting for previous observations, and then reconstructs the missing entries via least-squares estimation if required. Numerical examples are provided for direction-of-arrival estimation and matrix completion, comparing PETRELS with state of the art batch algorithms.


I. INTRODUCTION
Many real world datasets exhibit an embedding of lowdimensional structures in a high-dimensional manifold. When the embedding is assumed linear, the underlying lowdimensional structure becomes a linear subspace. Subspace Identification and Tracking (SIT) plays an important role in various signal processing tasks such as online identification of network anomalies [1], moving target localization [2], beamforming [3], and denoising [4]. Conventional SIT algorithms collect full measurements of the data stream at each time, and subsequently update the subspace estimate by utilizing Y. Chi  the track record of the stream history in different ways [5], [6].
Recent advances in Compressive Sensing (CS) [7], [8] and Matrix Completion (MC) [9], [10] have made it possible to infer data structure from highly incomplete observations. Compared with CS, which allows reconstruction of a single vector from only a few attributes by assuming it is sparse in a pre-determined basis or dictionary, MC allows reconstruction of a matrix from a few entries by assuming it is low rank. A popular method to perform MC is to minimize the nuclear norm of the corresponding matrix [9], [10] that the observed entries are satisfied. This method requires no prior knowledge of rank, in a similar spirit with ℓ 1 minimization [11] for sparse recovery in CS. Other approaches including greedy algorithms such as OptSpace [12] and ADMiRA [13] require an estimate of the matrix rank for initialization. Identifying the underlying low-rank structure in MC is equivalent to subspace identification in a batch setting. When the number of observed entries is slightly larger than the subspace rank, it has been shown that with high probability, it is possible to test whether a highly incomplete vector of interest lies in a known subspace [14]. Recent works on covariance matrix and principal components analysis of a dataset with missing entries also validate that it is possible to infer the principal components with high probability [15], [16].
In high-dimensional problems, it might be expensive and even impossible to collect data from all dimensions. For example in wireless sensor networks, collecting data from all sensors continuously will quickly drain the battery power. Ideally, we would prefer to obtain data from a fixed budget of sensors of each time to increase the overall battery life, and still be able to identify the underlying structure. Another example is in online recommendation systems, where it is impossible to expect rating feedbacks from all users on every product are available. Therefore it is of growing interest to identify and track a low-dimensional subspace from highly incomplete information of a data stream in an online fashion. In this setting, the estimate of the subspace is updated and tracked across time when new observations become available with low computational cost. The GROUSE algorithm [17] has been recently proposed for SIT from online partial observations using rank-one updates of the estimated subspace on the Grassmannian manifold. However, performance is limited by the existence of "barriers" in the search path [18] which result in GROUSE being trapped at a local minima. We demonstrate this behavior through numerical examples in Section VI in the context of direction-of-arrival estimation.
In this paper we further study the problem of SIT given partial observations from a data stream as in GROUSE. Our proposed algorithm is dubbed Parallel Estimation and Tracking by REcursive Least Squares (PETRELS). The underlying low-dimensional subspace is identified by minimizing the geometrically discounted sum of projection residuals on the observed entries per time index, via a recursive procedure with discounting for each row of the subspace matrix in parallel. The missing entries are then reconstructed via leastsquares estimation if required. The discount factor balances the algorithm's ability to capture long term behavior and changes to that behavior to improve adaptivity. We also benefit from the fact that our optimization of the estimated subspace is on all the possible low-rank subspaces, not restricted to the Grassmannian manifold. In the partial observation scenario, PETRELS always converges locally to a stationary point since it it a second-order stochastic gradient descent algorithm. In the full observation scenario, we prove that PETRELS actually converges to the global optimum by revealing its connection with the well-known Projection Approximation Subspace Tracking (PAST) algorithm [5]. Finally, we provide numerical examples to measure the impact of the discount factor, estimated rank and number of observed entries. In the context of direction-of-arrival estimation we demonstrate superior performance of PETRELS over GROUSE in terms of separating close-located modes and tracking changes in the scene. We also compare PETRELS with state of the art batch MC algorithms, showing it as a competitive alternative when the subspace is fixed.
The rest of the paper is organized as follows. Section II states the problem and provides background in the context of matrix completion and conventional subspace tracking. Section III describes the algorithm in details. Two extensions of PETRELS to improve robustness and reduce complexity are presented in Section IV. We discuss convergence issues of PETRELS in the full observation scenario in Section V. Section VI shows the numerical results and we conclude the paper in Section VII.

A. Problem Statement
We consider the following problem. At each time t, a vector x t ∈ R M is generated as: where the columns of U t ∈ R M×rt span a low-dimensional subspace, the vector a t ∈ R rt specifies the linear combination of columns and is Gaussian distributed as a t ∼ N (0, I rt ), and n t is an additive white Gaussian noise distributed as n t ∼ N (0, σ 2 I M ). The rank of the underlying subspace r t is not assumed known exactly and can be slowly changing over time. The entries in the vectors x t can be considered as measurements from different sensors in a sensor network, values of different pixels from a video frame, or movie ratings from each user.
We assume only partial entries of the full vector x t are observed, given by where ⊙ denotes point-wise multiplication, P t = diag{p t }, p t = [p 1t , p 2t , · · · , p Mt ] T ∈ {0, 1} M with p mt = 1 if the mth entry is observed at time t. We denote Ω t = {m : p mt = 1} as the set of observed entries at time t. In a random observation model, we assume the measurements are taken uniformly at random. We are interested in an online estimate of a low-rank subspace D n ∈ R M×r at each time index n, which identifies and tracks the changes in the underlying subspace, from streaming partial observations (y t , p t ) n t=1 . The rank of the estimated subspace D n is assumed known and fixed throughout the algorithm as r. In practice, we assume the upper bound of the rank of the underlying subspace U t is known, and sup t r t ≤ r. The desired properties for the algorithm include: • Low complexity: each step of the online algorithm at time index n should be adaptive with small complexity compared to running a batch algorithm using history data; • Small storage: The online algorithm should require a storage size that does not grow with the data size; • Convergence: The subspace sequence generated by the online algorithm should converge to the true subspace U t = U if it is constant. • Adaptivity: The online algorithm should be able to track the changes of the underlying subspace in a timely fashion.

B. Conventional Subspace Identification and Tracking
When x t 's are fully observed, our problem is equivalent to the classical SIT problem, which is widely studied and has a rich literature in the signal processing community. Here we describe the Projection Approximation Subspace Tracking (PAST) algorithm in details which is the closest to our proposed algorithm in the conventional scenario.
First, consider optimizing the scalar function with respect to a subspace W ∈ R M×r , given by be the data covariance matrix. It is shown in [5] that the global minima of (3) is the only stable stationary point, and is given by W = U r Q, where U r is composed of the r dominant eigenvectors of C x , and Q ∈ C r×r is a unitary matrix. Without loss of generality, we can choose U r = U. This motivates PAST to optimize the following function at time n without constraining W to have orthogonal columns: ≈ argmin where the expectation in (3) is replaced by geometrically reweighting the previous observations by α in (4), and is further approximated by replacing the second W by its previous estimate in (5). Based on (5), the subspace W n can be found by first estimating the coefficient vector a n using the previous subspace estimate as a n = W T n−1 x n , and updating the subspace as Suppose that α = 1 and denote R n = n i=1 a n a T n . In [19], the asymptotic dynamics of the PAST algorithm is described by its equilibrium as time goes to infinity using the Ordinary Differential Equation (ODE) below: whereã n = W T x n , R = R(t) and W = W(t) are continuous time versions of R n and W n , and † denotes the pseudo-inverse. It is proved in [19] that as t increases, W(t) converges to the global optima, i.e. to a matrix which spans the eigenvectors of C x corresponding to the r largest eigenvalues. In Section V we show that our proposed PETRELS algorithm becomes essentially equivalent to PAST when all entries of the data stream are observed, and can be shown to converge globally.
The PAST algorithm belongs to the class of power-based techniques, which include the Oja's method [20], the Novel Information Criterion (NIC) method [21] and etc: These algorithms are treated under a unified framework in [22] with slight variations for each algorithm. The readers are referred to [22] for details. In general, the estimate of the low-rank subspace W n ∈ R M×r is updated at time n as where C n is the sample data covariance matrix updated from and α n is a parameter between 0 and 1. The normalization in (7) assures that the updated subspace W n is orthogonal but this normalization is not performed strictly in different algorithms.
It is shown in [22] that these power-based methods guarantee global convergence to the principal subspace spanned by eigenvectors corresponding to the r largest eigenvalues of C x . If the entries of the data vector x t 's are fully observed, then C n converges to C x very fast, and this is exactly why the power-based methods perform very well in practice. When the data is highly incomplete, the convergence of (8) is very slow since only a small fraction |Ω n | 2 /n 2 of entries in C n−1 are updated, where |Ω n | is the number of observed entries at time n, making direct adoption of the above method unrealistic in the partial observation scenario.

C. Matrix Completion
When only partial observations are available and U t = U are fixed, our problem is closely related to the Matrix Completion (MC) problem, which has been extensively studied recently. Assume X ∈ R M×n is a low-rank matrix, P is a binary M × n mask matrix with 0 at missing entries and 1 at observed entries. Let Y = P ⊙ X = [y 1 , . . . , y n ] be the observed partial matrix where the missing entries are filled in as zero, and ⊙ denotes point-wise multiplication. MC aims to solve the following problem: i.e. to find a matrix with the minimal rank such that the observed entries are satisfied. This problem is combinatorially intractable due to the rank constraint. It has been shown in [9] that by replacing the rank constraint with nuclear norm minimization, (9) can be solved by a convex optimization problem, resulting in the following spectralregularized MC problem: where Z * is the nuclear norm of Z, i.e. the sum of singular values of Z, and µ > 0 is a regularization parameter. Under mild conditions, the solution of (10) is the same as that of (9) [9]. The nuclear norm [23] of Z is given by where U ∈ C M×r and V ∈ C n×r . Substituting (11) in (10) we can rewrite the MC problem as Our problem formulation can be viewed as an online way of solving the above batch-setting MC problem. Consider a random process {n t } where each n t is drawn uniformly from {1, . . . , n}, and a data stream is constructed where the data at each time is given as x nt , i.e. the n t th column of X. Compared with (1), the subspace is fixed as U t = U since we draw columns from a fixed low-rank matrix. Each time we only observe partial entries of x nt , given as y nt = p nt ⊙x nt , where P nt is the n t th column of P. The problem of MC becomes equivalent to retrieving the underlying subspace U from the data stream (y nt , p nt ) ∞ t=1 . After estimating U, the low-rank matrix X can be recovered via least-squares estimation. The online treatment of the batch MC problem has potential advantages for avoiding large matrix manipulations. We will compare the PETRELS algorithm against some of the popular MC methods in Section VI.

III. THE PETRELS ALGORITHM
We now describe our proposed Parallel Estimation and Tracking by REcursive Least Squares (PETRELS) algorithm.

A. Objective Function
We first define the function f t (D) at each time t = 1, · · · , n for a fixed subspace D ∈ R M×r , which is the total projection residual on the observed entries, Here r is the rank of the estimated subspace, which is assumed known and fixed throughout the algorithm 1 . We aim to minimize the following loss function at each time n with respect to the underlying subspace: where D n is the estimated subspace of rank r at time n, and the parameter 0 ≪ λ ≤ 1 discounts past observations. Before developing PETRELS we note that if there are further constraints on the coefficients a t 's, a regularization term can be incorporated as: where p ≥ 0. For example, p = 1 enforces a sparse constraint on a t , and p = 2 enforces a norm constraint on a t . In (14) the discount factor λ is fixed, and the influence of past estimates decreases geometrically; a more general online objective function can be given as where the sequence {λ n } is used to control the memory and adaptivity of the system in a more flexible way.
To motivate the loss function in (14) we note that if U t = U is not changing over time, then the RHS of (14) is minimized to zero when D n spans the subspace defined by U. If U t is slowly changing, then λ is used to control the memory of the system and maintain tracking ability at time n. For example, by using λ → 1 the algorithm gradually loses its ability to forget the past.
Fixing D, f t (D) can be written as Plugging this back to (14) the exact optimization problem becomes: This problem is difficult to solve over D and requires storing all previous observations. Instead, we propose PETRELS to approximately solve this optimization problem.

B. PETRELS
The proposed PETRELS algorithm, as summarized by Algorithm 1, alternates between coefficient estimation and subspace update at each time n. In particular, the coefficient vector a n is estimated by minimizing the projection residual on the previous subspace estimate D n−1 : a n = argmin 1 The rank may not equal the true subspace dimension.

Algorithm 1 PETRELS for SIT from Partial Observations
Input: a stream of vectors y t and observed pattern P t . Initialization: an M ×r random matrix D 0 , and (R 0 m ) † = δI r , δ > 0 for all m = 1, · · · , M . 1: for n = 1, 2, · · · do 2: a n = (D T n−1 P n D n−1 ) † D T n−1 y n . 3:x n = D n−1 a n . 4: for m = 1, · · · , M do 5: where D 0 is a random subspace initialization. The full vector x n is then estimated as: The subspace D n is then updated by minimizing where a t , t = 1, · · · , n are estimates from (18). Comparing (20) with (14), the optimal coefficients are substituted for the previous estimated coefficients. This results in a simpler problem for finding D n . The discount factor mitigates the error propagation and compensates for the fact that we used the previous coefficients updated rather than solving (14) directly, therefore improving the performance of the algorithm. The objective function in (20) can be equivalently decomposed into a set of smaller problems for each row of for m = 1, · · · , M . To find the optimal d n m , we equate the derivative of (21) to zero, resulting in

This equation can be rewritten as
where R n m = n t=1 λ n−t p mt a t a T t and s n m = n t=1 λ n−t p mt x mt a t . Therefore, d n m can be found as When R n m is not invertible, (23) is the least-norm solution to d n m . We now show how (22) can be updated recursively. First we rewrite R n m = λR n−1 m + p mn a n a T n , for all m = 1, · · · , M . Then we plug (24) and (25) into (22), and get R n m d n m = λs n−1 m + p mn x mn a n = λR n−1 m d n−1 m + p mn x mn a n = R n m d n−1 m − p mn a n a T n d n−1 m + p mn x mn a n = R n m d n−1 where d n−1 m is the row estimate in the previous time n − 1. This results in a parallel procedure to update all rows of the subspace matrix D n , give as Finally, by the Recursive Least-Squares (RLS) updating formula for the general pseudo-inverse matrix [24], [25], (R n m ) † can be easily updated without matrix inversion using Here G n m = (β n m ) −1 v n m (v n m ) T , with β n m and v n m given as To enable the RLS procedure, the matrix (R 0 m ) † is initialized as a matrix with large entries on the diagonal, which we choose arbitrarily as the identity matrix (R 0 m ) † = δI r , δ > 0 for all m = 1, · · · , M . It is worth noting that implementation of the fast RLS update rules is in general very efficient. However, caution needs to be taken since direct application of fast RLS algorithms suffer from numerical instability of finite-precision operations when running for a long time [26].

C. Second-Order Stochastic Gradient Descent
The PETRELS algorithm can be regarded as a second-order stochastic gradient descent method to solve (14) by using d n−1 m , m = 1, · · · , M as a warm start at time n. Specifically, we can write the gradient of f n (D) in (13) at D n−1 as where a n is given in (18). Then the gradient of F n (D) at D n−1 is given as The Hessian for each row of D at d n−1 m is therefore It follows that the update rule for each row d m given in (27) can be written as which is equivalent to second-order stochastic gradient descent. Therefore, PETRELS converges to a stationary point of F n (D) [27], [28]. Compared with first-order algorithms, PETRELS enjoys a faster convergence speed to the stationary point [27], [28].

D. Comparison with GROUSE
The GROUSE algorithm [17] proposed by Balzano et. al. addresses the same problem of online identification of low-rank subspace from highly incomplete information. The GROUSE method can be viewed as optimizing (14) for λ = 1 at each time n using a first-order stochastic gradient descent on the orthogonal Grassmannian defined as G r = {D ∈ R M×r : D T D = I r } instead of R M×r . Thus, GROUSE aims to solve the following optimization problem, GROUSE updates the subspace estimate along the direction of ∇f t (D)| D=Dn−1 on G r , given by where σ = x t 2 r t 2 , and η n is the step-size at time n. At each step GROUSE also alternates between coefficient estimation (18) and subspace update (33). Moreover, the resulting algorithm is a fast rank-one update on D n−1 at each time n. Given that it is a first-order gradient descent algorithm, convergence to a stationary point but not global optimal is guaranteed under mild conditions on the step-size. Specifically, if the step size satisfies then GROUSE is guaranteed to converge to a stationary point of G n (D). However, due to the existence of "barriers" in the search path on the Grassmannian [18], GROUSE may be trapped at a local minima as shown in Section VI in the example of direction-of-arrival estimation. Although both PETRELS and GROUSE have a tuning parameter, compared with the step-size in GROUSE, the discount factor in PE-TRELS is an easier parameter to tune. For example, without discounting (i.e. λ = 1) PETRELS can still converge to the global optimal given full observations as shown in Section V, while this is impossible to achieve with a first-order algorithm like GROUSE if the step size is not tuned properly to satisfy (34).
If we relax the objective function of GROUSE (32) to all rank-r subspaces R M×r , given as then the objective function becomes equivalent to PETRELS without discounting. It is possible to use a different formula-tion of second-order stochastic gradient descent with step-size to solve (35), yielding the update rule for each row of D n as where H n (d n−1 m , λ = 1) is given in (30), and γ n is the stepsize at time n. Compared with the update rule for PETRELS in (31), the discount parameter has a similar role as the step-size, but weights the contribution of previous data input geometrically. However, in this paper we didn't investigate the performance of this alternative update rule in (36).

E. Complexity Issues
We compare both storage complexity and computational complexity for PETRELS, GROUSE and the PAST algorithm. The storage complexity of PAST and GROUSE is O(M r), which is the size of the low-rank subspace. On the other hand, PETRELS has a larger storage complexity of O(M r 2 ), which is the total size of R n m 's for each row. In terms of computational complexity, PAST has a complexity of O(M r), while PETRELS and GROUSE have a similar complexity on the order of O(|Ω t |r 2 ), where the main complexity comes from computation of the coefficient (18). This indicates another merit of dealing with partial observations, i.e. to reduce computational complexity when the dimension is high.

A. Simplified PETRELS
In the subspace update step of PETRELS in (20), consider replacing the objective function in (14) by where a t andx t , t = 1, · · · , n are estimates from earlier steps in (18) and (19). The only change we made is to remove the partial observation operator from the objective function, and replace it by the full vector estimate. It remains true that d n m = argmin dmFn (d m ) = d n−1 m if the corresponding mth entry of x n is unobserved, i.e. m / ∈ Ω n , sincê ∈ Ω n . This modification leads to a simplified update rule for R n m , since now the updating formula for all rows d m 's is the same, where R n m = R n = λR n−1 + a n a T n for all m. The row updating formula (27) is replaced by D n = D n−1 + P n (x n − D n−1 a n )a T n R † n , which further saves storage requirement for the PETRELS algorithm from O(M r 2 ), to O(M r) which is the size of the subspace. We compared the performance of the simplified PETRELS against PETRELS in Section VI, which converges slower than PETRELS but might have an advantage if the subspace rank is underestimated.

B. Incorporating Prior Information
It is possible to incorporate regularization terms into PE-TRELS to encode prior information about the data stream.
Here we outline the regularization on the subspace D in the subspace update step, such that at each time n, D n is updated via where µ n > 0 is the regularization parameter. Similar as PETRELS, (39) can be decomposed for each row of D = The matrix T n m can be updated as T n m = λT n−1 m + p mn a t a T t + (µ n − λµ n−1 )I r , and s n m can be updated as (25). However the fast RLS algorithm no longer applies here, so additional complexity for matrix inversion is required. It is worth noticing that (39) closely resembles the matrix completion formula (12) when V is fixed and composed of columns of a t 's.

V. GLOBAL CONVERGENCE WITH FULL OBSERVATION
In the partial observation regime, the PETRELS algorithm always converges to a stationary point of F n (D), given it's a second-order stochastic gradient descent method in Section III-C, but whether it converges to the global optimal remains open. However, in the full observation regime, i.e. y n = x n for all n, we can show that the PETRELS algorithm converge globally as below.
In this case, PETRELS becomes essentially equivalent to the conventional PAST algorithm [5] for SIT except that the coefficient is estimated differently. Specifically, in PAST it is estimated as a n = D T n−1 y n = D T n−1 x n , while in PETRELS it is estimated as a n = (D T n−1 D n−1 ) −1 D T n−1 x n . Now let λ = 1, similar to PAST in [19], the asymptotic dynamics of the PETRELS algorithm can be described by the ODE below, for any function of D T D. Hence, Therefore (40) and (41) can be rewritten aṡ which is equivalent to the ODE of PAST. Hence we conclude that PETRELS will converge to the global optima in the same dynamic as the PAST algorithm.

VI. NUMERICAL RESULTS
Our numerical results fall into four parts. First we examine the influence of parameters specified in the PETRELS algorithm, such as discount factor, rank estimation, and its robustness to noise level. Next we look at the problem of direction-of-arrival estimation and show PETRELS demonstrates performance superior to GROUSE by identifying and tracking all the targets almost perfectly even in low SNR. Thirdly, we compare our approach with matrix completion, and show that PETRELS is at least competitive with state of the art batch algorithms. Finally, we provide numerical simulations for the extensions of the PETRELS algorithm.

A. Choice of Parameters
At each time t, a vector x t is generated as where D true is an r-dimensional subspace generated with i.i.d. N (0, 1) entries, a t is an r×1 vector with i.i.d. N (0, 1) entries, and n t is an m × 1 Gaussian noise vector with i.i.d. N (0, ǫ 2 ) entries. We further fix the signal dimension m = 500 and the subspace rank r true = 10. We assume that a fixed number of entries in x t , denoted by K, are revealed each time. This restriction is not necessary for the algorithm to work as shown in matrix completion simulations, but we make it here in order to get a meaningful estimate of a t . Denoting the estimated subspace byD, we use the normalized subspace reconstruction error to examine the algorithm performance. This is calculated where PD ⊥ is the projection operator to the orthogonal subspaceD ⊥ .
The choice of discount factor λ plays an important role in how fast the algorithm converges. With K = 50, a mere 10% percent of the full dimension, the rank is given accurately as r = 10 in a noise-free setting where ǫ = 0. We run the algorithm to time n = 2000 for the same data, and find that the normalized subspace reconstruction error is minimized when λ is around 0.98 as shown in Fig. 1. Hence, we will keep λ = 0.98 hereafter.  In reality it is almost impossible to accurately estimate the intrinsic rank in advance. Fortunately the convergence rate of our algorithm degrades gracefully as the rank estimation error increases. In Fig. 2, the evolution of normalized subspace error is plotted against data stream index, for rank estimation r = 10, 12, 14, 16, 18. We only examine over-estimation of the rank here since this is usually the case in applications.
In the next section we show examples for the case of rank underestimation.
Taking more measurements per time leads to faster convergence since it is approaching the full information regime, as shown in Fig. 3. Theoretically it requires M ∼ O(r log r) ≈ 23 measurements to test if an incomplete vector is within a subspace of rank r [14]. The simulation shows our algorithm can work even when M is close to this lower bound.
Finally the robustness of the algorithm is tested against the noise variance ǫ 2 in Fig. 4, where the normalized subspace error is plotted against data stream index for different noise levels ǫ. The estimated subspace deviates from the ground truth as we increase the noise level, hence the normalized subspace error degrades gracefully and converges to an error floor determined by the noise variance.
We now consider a scenario where a subspace of rank r = 10 changes abruptly at time index n = 3000 and n = 5000, and examine the performance of GROUSE [17] and PETRELS in Fig. 5 when the rank is over-estimated by 4 and the noise level is ǫ = 10 −3 . The normalized residual error for data stream, calculated as P n (x n −x n ) 2 / P n x n 2 , is shown in Fig. 5 (a), and the normalized subspace error is shown in Fig. 5 (b) respectively. Both PETRELS and GROUSE can successfully track the changed subspace, but PETRELS can track the change faster.

B. Direction-Of-Arrival Analysis
Given GROUSE [17] as a baseline, we evaluate the resilience of our algorithm to different data models and applications. We use the following example of Direction-Of-Arrival analysis in array processing to compare the performance of these two methods. Assume there are n = 256 sensors from a linear array, and the measurements from all sensors at time t are given as (43) Here V ∈ C n×p is a Vandermonde matrix given by where α i (ω i ) = [1, e j2πωi , · · · , e j2πωi(n−1)] ] T , 0 ≤ ω i < 1; Σ = diag{d} = diag{d 1 , · · · , d p } is a diagonal matrix which characterizes the amplitudes of each mode. The coefficients a t are generated with N (0, 1) entries, and the noise is generated with N (0, ǫ 2 ) entries, where ǫ = 0.1. Each time we collect measurements from K = 30 random sensors. We are interested in identifying all . This can be done by applying the well-known ESPRIT algorithm [29] to the estimated subspaceD of rank r, where r is specified a-priori corresponding to the number of modes to be estimated. Specifically, if D 1 =D(1 : n − 1) and D 2 =D(2 : n) are the first and the last n − 1 rows ofD, then from the eigenvalues of the matrix T = D † 1 D 2 , denoted by λ i , i = 1, · · · , r, the set of {ω i } p i=1 can be recovered as The ESPRIT algorithm also plays a role in recovery of multipath delays from low-rate samples of the channel output [30]. We show that in a dynamic setting when the underlying subspace is varying, PETRELS does a better job of discarding out-of-date modes and picking up new ones in comparison with GROUSE. We divide the running time into 4 parts, and the frequencies and amplitudes are specified as follows:   6 shows the ground truth of mode locations and amplitudes for the scenario above. Note that there are three closely located modes and one weak mode in the beginning, which makes the task challenging. We compare the performance of PETRELS and GROUSE. The rank specified in both algorithms is r = 10, which is the number of estimated modes at each time index; in our case it is twice the number of true modes. 2 Each time both algorithms estimated 10 modes, with their amplitude shown shown against the data stream index in Fig. 7  (a) and (b). The color shows the amplitude corresponding to the color bar. The direction-of-arrival estimations in Fig. 7 (a) and (b) are further thresholded with respect to level 0.5, and the thresholded results are shown in Fig. 7 (c) and (d) for PETRELS and GROUSE respectively. PETRELS identifies all modes correctly. In particular PETRELS distinguishes the three closely-spaced modes perfectly in the beginning, and identifies the weak modes that come in the scene at a later time. With GROUSE the closely spaced nodes are erroneously estimated as one mode, the weak mode is missing, and spurious modes have been introduced. PETRELS also fully tracked the later changes in accordance with the entrance and exit of each mode, while GROUSE is not able to react to changes in the data model.
Since the number of estimated modes at each time is greater than the number of true modes, the additional rank in the estimated subspace contributes "auxiliary modes" that do not belong to the data model. In PETRELS these modes exhibit as scatter points with small amplitudes as in Fig. 7 (a), so they will not be identified as actual targets in the scene. While in GROUSE these auxiliary modes are tracked and appear as spurious modes. All changes are identified and tracked successfully by PETRELS, but not by GROUSE.

C. Matrix Completion
We next compare performance of PETRELS for matrix completion against batch algorithms LMaFit [32], FPCA [33], Singular Value Thresholding (SVT) [34], OptSpace [12] and GROUSE [17]. The low-rank matrix is generated from a matrix factorization model with X = UV T ∈ R 1000×2000 , where U ∈ R 1000×10 and V ∈ R 2000×10 , all entries in U and V are generated from standard normal distribution N (0, 1) (Gaussian data) or uniform distribution U[0, 1] (uniform data). The sampling rate is taken to be 0.05, so only 5% of all entries are revealed. The running time is plotted against the normalized matrix reconstruction error, calculated as X − X F / X F , wherê X is the reconstructed low-rank matrix for Gaussian data and uniform data respectively in Fig. 8 (a) and (b). PETRELS matches the performance of batch algorithms on Gaussian data and improves upon the accuracy of most algorithms on uniform data, where the Grassmaniann-based optimization approach may encounter "barriers" for its convergence. Note that different algorithms have different input parameter requirements. For example, OptSpace needs to specify the tolerance to terminate the iterations, which directly decides the trade-off between accuracy and running time; PETRELS and GROUSE require an initial estimate of the rank. Our simulation here only shows one particular realization and we simply conclude that PETRELS is competitive.

D. Simplified PETRELS
Under the same simulation setup as for Fig. 2 except that the subspace of rank 10 is generated byD true = D true Σ, where Σ is a diagonal matrix with 5 entries from N (0, 1) and 5 entries from 0.01 · N (0, 1), we examine the performance of the simplified PETRELS algorithm (with optimized λ = 0.9) in Section IV A and the original PETRELS (with λ = 0.98) algorithm when the rank of the subspace is over-estimated as 12 or under-estimated as 8. When the rank of D is over-estimated, the change in (9) will introduce more errors and converges slower compared with the original PETRELS algorithm; however, when the rank of D is under-estimated, the simplified PETRELS performs better than PETRELS. This is an interesting feature of the proposed simplification, and quantitative justification of this phenomenon is beyond  the scope of this paper. Intuitively, when the rank is underestimated, the simplified PETRELS also uses the interpolated entries to update the subspace estimate, which seems to help the performance.

VII. CONCLUSIONS
We considered the problem of reconstructing a data stream from a small subset of its entries, where the data stream is assumed to lie in a low-dimensional linear subspace, possibly corrupted by noise. This has significant implications for lessening the storage burden and reducing complexity, as well as tracking the changes for applications such as video denoising, network monitoring and anomaly detection when the problem size is large. The well-known low-rank matrix completion problem can be viewed as a batch version of our problem. The PETRELS algorithm first identifies the underlying lowdimensional subspace via a discounted recursive procedure for each row of the subspace matrix in parallel, then reconstructs the missing entries via least-squares estimation if required. The discount factor allows the algorithm to capture long-term behavior as well as track the changes of the data stream. We show that PETRELS converges to a stationary point given it is a second-order stochastic gradient descent algorithm. In the full observation scenario, we further prove that PETRELS actually convergence globally by revealing its connection with the PAST algorithm. We demonstrate superior performance of PETRELS in direction-of-arrival estimation and showed that it is competitive with state of the art batch matrix completion algorithms.