Provable Dynamic Robust PCA or Robust Subspace Tracking

TL;DR

This paper introduces simple-ReProCS, a novel online algorithm for dynamic robust PCA that effectively tracks changing low-dimensional subspaces in data with sparse outliers, under weaker assumptions than previous methods.

Contribution

It provides the first theoretical guarantee for dynamic RPCA under relaxed assumptions, improving outlier tolerance and demonstrating online, fast, and memory-efficient performance.

Findings

Guarantees dynamic RPCA under weakened assumptions.

Achieves higher outlier tolerance than previous methods.

Operates online with near-optimal memory complexity.

Abstract

Dynamic robust PCA refers to the dynamic (time-varying) extension of robust PCA (RPCA). It assumes that the true (uncorrupted) data lies in a low-dimensional subspace that can change with time, albeit slowly. The goal is to track this changing subspace over time in the presence of sparse outliers. We develop and study a novel algorithm, that we call simple-ReProCS, based on the recently introduced Recursive Projected Compressive Sensing (ReProCS) framework. Our work provides the first guarantee for dynamic RPCA that holds under weakened versions of standard RPCA assumptions, slow subspace change and a lower bound assumption on most outlier magnitudes. Our result is significant because (i) it removes the strong assumptions needed by the two previous complete guarantees for ReProCS-based algorithms; (ii) it shows that it is possible to achieve significantly improved outlier tolerance, compared with all existing RPCA or dynamic RPCA solutions by exploiting the above two simple extra assumptions; and (iii) it proves that simple-ReProCS is online (after initialization), fast, and, has near-optimal memory complexity.

Figures

Tables

Table 1. TABLE I: Comparing s-ReProCS with other RPCA solutions with complete guarantees. For simplicity, we ignore all dependence on condition numbers. In this table r L subscript 𝑟 𝐿 r_{\scriptscriptstyle{L}} is the rank of the entire matrix 𝑳 𝑳 \bm{L} , while r 𝑟 r is the maximum rank of any sub-matrices of consecutive columns of 𝑳 𝑳 \bm{L} of the form 𝑳 [ t j , t j + 1 ) subscript 𝑳 subscript 𝑡 𝑗 subscript 𝑡 𝑗 1 \bm{L}_{[t_{j},t_{j+1})} and thus r ≤ r L 𝑟 subscript 𝑟 𝐿 r\leq r_{\scriptscriptstyle{L}} . We show the unrealistic assumptions in red .

Algorithm	Outlier tolerance	Assumptions	Memory, Time,	# params.
PCP (C) [2]	$\text{max-outlier-frac-row}\in 𝒪\left(1\right)$	outlier support: uniform random	Memory: $𝒪\left(nd\right)$	zero
(offline)	$\text{max-outlier-frac-col}\in 𝒪\left(1\right)$	$r_L\le c\min (n,d)/{\log }^{2}n$	Time: $𝒪\left(n{d}^2\frac{1}{ϵ}\right)$
PCP (H) [4]	$\text{max-outlier-frac-row}\in 𝒪\left(1/r_L\right)$		Memory: $𝒪\left(nd\right)$	$2$
(offline)	$\text{max-outlier-frac-col}\in 𝒪\left(1/r_L\right)$		Time: $𝒪\left(n{d}^2\frac{1}{ϵ}\right)$
AltProj [6],	$\text{max-outlier-frac-row}=O\left(1/r_L\right)$		Memory: $𝒪\left(nd\right)$	2
(offline)	$\text{max-outlier-frac-col}\in O\left(1/r_L\right)$		Time: $𝒪\left(nd{r}_L^2\log \frac{1}{ϵ}\right)$
RPCA-GD [8]	$\text{max-outlier-frac-row}\in 𝒪\left(1/r_L^{1.5}\right)$		Memory: $𝒪\left(nd\right)$	5
(offline)	$\text{max-outlier-frac-col}\in 𝒪\left(1/r_L^{1.5}\right)$		Time: $𝒪\left(nd{r}_L\log \frac{1}{ϵ}\right)$
NO-RMC [18]	$\text{max-outlier-frac-row}\in O\left(1/r_L\right)$	$Cn\ge d\ge cn$	Memory: $𝒪\left(nd\right)$	3
(offline)	$\text{max-outlier-frac-col}\in 𝒪\left(1/r_L\right)$		Time: $𝒪\left(n{r}_L^3{\log }^{2}n{\log }^2\frac{1}{ϵ}\right)$
s-ReProCS	${\text{max-outlier-frac-row}}^{𝜶}\in 𝓞\left(\mathrm{𝟏}\right)$	most outlier magnitudes lower bounded	Memory: $𝒪\left(nr\log n\right)$	4
(online)	$\text{max-outlier-frac-col}\in 𝓞\left(\mathrm{𝟏}/𝒓\right)$	slow subspace change	Time: $𝒪\left(ndr\log \frac{1}{ϵ}\right)$
(this work)		first $Cr$ samples: AltProj assumptions	Detect delay: $2\alpha =Cr\log n$
			Tracking Delay: $K\alpha =Cr\log n\log \left(1/ϵ\right)$

Table 2. TABLE II: Comparing s-ReProCS with online or tracking approaches for RPCA. We show the unrealistic assumptions in red . Here, f 𝑓 f denotes the condition number of 𝚲 𝚲 \bm{\Lambda} , r 𝑟 r is the maximum dimension of the subspace at any time, and r L subscript 𝑟 𝐿 r_{\scriptscriptstyle{L}} refers to the rank of matrix 𝑳 𝑳 \bm{L} . Thus r ≤ r L 𝑟 subscript 𝑟 𝐿 r\leq r_{\scriptscriptstyle{L}} . Here, s-ReProCS-no-delete refers to Algorithm 4 without the subpace deletion step.

Algorithm	Outlier tolerance	Assumptions	Memory, Time
orig-ReProCS [7, 12]	${\text{max-outlier-frac-row}}^{\alpha }\in 𝒪\left(1/f^2\right)$	outlier support: moving object model,	Memory: $O\left(n{r}^2/{ϵ}^2\right)$
(online)	$\text{max-outlier-frac-col}\in O\left(1/r_L\right)$	unrealistic subspace change model,	Time: $𝒪\left(ndr\log \frac{1}{ϵ}\right)$
		changed eigenvalues small for some time,	Detect Delay: $2\alpha =\frac{C{r}^2\log n}{{ϵ}^2}$
		outlier mag. lower bounded,	Tracking Delay: $K\alpha =\frac{C{r}^2\log n\log \left(1/ϵ\right)}{{ϵ}^2}$
		$x_{\min }\ge 14\left[c{\gamma }_{\mathrm{new}}+\sqrt{\tilde{\epsilon }}\left(\sqrt{r}+\sqrt{c}\right)\right]$
		where, ${\gamma }_{\mathrm{new}}$ quantifies slow subspace change
		init data: AltProj assumptions,
		$d\ge C{r}^2/{ϵ}^2$
Modified-PCP [13]	${\text{max-outlier-frac-row}}^{\alpha }\in 𝒪\left(1\right)$	outlier support: uniform random	Memory: $𝒪\left(nr{\log }^{2}n\right)$
(piecewise batch)	$\text{max-outlier-frac-col}\in 𝒪\left(1\right)$	unrealistic subspace change model	Time: $𝒪\left(\frac{ndr{\log }^{2}n}{ϵ}\right)$
		$r_L\le c\min (n,d)/{\log }^{2}n$	Detect delay: $\infty$
ORPCA [14]	Has a partial guarantee – assumes algorithm estimates at each time $t$ are full rank
GRASTA [15]	Has no theoretical guarantees
s-ReProCS	${\text{max-outlier-frac-row}}^{𝜶}\in 𝓞\left(\mathrm{𝟏}/{𝒇}^{\mathrm{𝟐}}\right)$	most outlier magnitudes lower bounded	Memory: $𝒪\left(nr\log n\right)$
(online)	$\text{max-outlier-frac-col}\in 𝓞\left(\mathrm{𝟏}/𝒓\right)$	$x_{\min }\ge 15C\left(2\tilde{\epsilon }\sqrt{r{\lambda }^{+}}+\Delta \sqrt{{\lambda }_{\mathrm{ch}}}\right)$	Time: $𝒪\left(ndr\log \frac{1}{ϵ}\right)$
(this work)		slow subspace change	Detect delay: $2\alpha =Cr\log n$
		first $Cr$ samples: AltProj assumptions	Tracking Delay: $K\alpha =Cr\log n\log \left(1/ϵ\right)$
s-ReProCS-no-delete	${\text{max-outlier-frac-row}}^{𝜶}\in 𝓞\left(\mathrm{𝟏}\right)$	most outlier magnitudes lower bounded	Memory: $𝒪\left(nr\log n\right)$
(online)	$\text{max-outlier-frac-col}\in 𝓞\left(\mathrm{𝟏}/{𝒓}_{𝑳}\right)$	$x_{\min }\ge 15C\left(2\tilde{\epsilon }\sqrt{r{\lambda }^{+}}+\Delta \sqrt{{\lambda }_{\mathrm{ch}}}\right)$	Time: $𝒪\left(ndr\log \frac{1}{ϵ}\right)$
(this work)		slow subspace change	Detect delay: $2\alpha =Cr\log n$
		first $Cr$ samples: AltProj assumptions	Tracking Delay: $K\alpha =Cr\log n\log \left(1/ϵ\right)$
ReProCS-NORST	$\text{max-outlier-frac-row}=𝒪\left(1/f^2\right)$	outlier mag. lower bounded	Memory: $𝒪\left(nr\log n\log \frac{1}{ϵ}\right)$
[20, 11] (online)	$\text{max-outlier-frac-col}=𝒪\left(1/r\right)$	$x_{\min }\ge C_1\sqrt{r{\lambda }^{+}}\left(\Delta +2\tilde{\epsilon }\right)$	Time: $𝒪\left(ndr\log \frac{1}{ϵ}\right)$
(follow-up		slow subspace change or fixed subspace	Detect delay: $Cr\log n$
to this work)		first $Cr$ samples: AltProj assumptions	Tracking Delay: $Cr\log n\log \left(1/ϵ\right)$

Table 3. TABLE III: List of Symbols and Assumptions used in the Main Result 2.2 , and Corollary 2.3 . (Note: We show that whp, t ^ j ≥ t j subscript ^ 𝑡 𝑗 subscript 𝑡 𝑗 {\hat{t}}_{j}\geq t_{j} and t ^ j + ( K + 1 ) ​ α ≤ t j + 1 subscript ^ 𝑡 𝑗 𝐾 1 𝛼 subscript 𝑡 𝑗 1 {\hat{t}}_{j}+(K+1)\alpha\leq t_{j+1} and hence, whp, 𝒥 0 , 𝒥 K + 2 subscript 𝒥 0 subscript 𝒥 𝐾 2 \mathcal{J}_{0},\mathcal{J}_{K+2} are non-empty intervals.

Observations: ${𝒚}_t={\ell }_t+{𝒙}_t+{𝒗}_t$, where, for $t\in [t_j,t_{j+1})$, ${𝒯}_t$ is support of ${𝒙}_t$.
Subspace Change		Principal Subspace Coefficients, ${𝒂}_t$’s
${𝑷}_{\left(t\right)}={𝑷}_{\left(t_j\right)}={𝑷}_j\text{ for all }t\in [t_j,t_{j+1}),j=0,1,\dots ,J$		element-wise bounded, zero mean,
${𝑷}_{\mathrm{ch}}\equiv {𝑷}_{j-1,\mathrm{ch}}$	Changing direction from $\mathrm{span}\left({𝑷}_{j-1}\right)$ at $t_j$	mutually independent with identical covariance (See Sec. II-B)
${𝑷}_{\mathrm{fix}}\equiv {𝑷}_{j-1,\mathrm{fix}}$	Fixed directions from $\mathrm{span}\left({𝑷}_{j-1}\right)$ at $t_j$
${𝑷}_{\mathrm{new}}\equiv {𝑷}_{j-1,\mathrm{new}}$	New direction from $\mathrm{span}\left({𝑷}_{j-1,⟂}\right)$ added at $t_j$	${\lambda }^{+}$	${\lambda }_{\max }\left(𝚲\right)$
${𝑷}_{\mathrm{rot}}\equiv {𝑷}_{j-1,\mathrm{rot}}$	Rotated version of ${𝑷}_{\mathrm{ch}}$	${\lambda }^{-}$	${\lambda }_{\min }\left(𝚲\right)$
$\mathrm{SE}({𝑷}_{j-1},{𝑷}_j)=\mathrm{SE}({𝑷}_{j-1,\mathrm{ch}},{𝑷}_{j,\mathrm{rot}})\le \Delta$		$f:={\lambda }^{+}/{\lambda }^{-}$	Condition Number of $𝚲$
${𝑷}_{j,\mathrm{new}}:=\frac{(𝑰-{𝑷}_{j-1,\mathrm{ch}}{𝑷}_{j-1,\mathrm{ch}}{}_{}{}^{\prime }){𝑷}_{j,\mathrm{rot}}}{\mathrm{SE}({𝑷}_{j-1,\mathrm{ch}},{𝑷}_{j,\mathrm{rot}})}$
See (5) for equivalent generative model.
${\max }_j{\max }_i\parallel \mathrm{basis}{\left([{𝑷}_{j-1},{𝑷}_j]\right)}^i\parallel \le \sqrt{\frac{\mu \left(r+1\right)}{n}}$ which implies (13) holds.
Outliers		Intervals for $j$-th subspace change and tracking
$x_{\min ,t}:={\min }_{i\in {𝒯}_t}\left|{\left({𝒙}_t\right)}_i\right|$	Min. outlier magnitude at $t$	${𝒥}_0:=[t_j,{\hat{t}}_j)$	interval before change detected
$x_{\min }:={\min }_t{x}_{\min ,t}$	Min. outlier magnitude	${𝒥}_k:=[{\hat{t}}_j+\left(k-1\right)\alpha ,{\hat{t}}_j+k\alpha )$	$k$-th subspace update interval
$s:=\text{max-outlier-frac-col}\cdot n$	Cardinality of support set of ${𝒙}_t$	${𝒥}_{K+1}:=[{\hat{t}}_j+K\alpha ,{\hat{t}}_j+\left(K+1\right)\alpha )$	SVD-re-estimation interval
${\text{max-outlier-frac-row}}^{\alpha }\le {\rho }_{\mathrm{row}}$	See (8)	${𝒥}_{K+2}:=[{\hat{t}}_j+\left(K+1\right)\alpha ,t_{j+1})$	Final interval
${\rho }_{\mathrm{row}}=0.01/f^2$
$\text{max-outlier-frac-col}\le {\rho }_{\mathrm{col}}=\frac{0.01}{2\mu \left(r+1\right)}$	See Theorem 2.2

Table 4. TABLE IV: List of symbols and their associated meaning for understanding the proof of Theorems 2.2 and 5.7 . The complete definitions can be found in Definitions 5.12 and A.20 . We also provide the location of the proof for each of events/scalars where applicable in parenthesis.

Symbol	Meaning
Preliminaries (Stated in Definitions 5.12 and A.20)
$\theta :={\theta }_j$	Angle of $j$-th subspace change
${𝑷}_{\ast }:={𝑷}_{j-1}$, ${𝑷}_{\mathrm{new}}:={𝑷}_{j,\mathrm{new}}$, ${𝑷}_{\mathrm{ch}}:={𝑷}_{j-1,\mathrm{ch}}$, ${𝑷}_{\mathrm{fix}}:={𝑷}_{j-1,\mathrm{fix}}$	Parts of the $j$-th subspace
${𝑷}_{\mathrm{rot}}:={𝑷}_{j,\mathrm{rot}}:=\left({𝑷}_{\mathrm{ch}}\cos \theta +{𝑷}_{\mathrm{new}}\sin \theta \right)$, $𝑷:={𝑷}_j$
${\widehat{𝑷}}_{\ast }:={\widehat{𝑷}}_{j-1}$, $\widehat{𝑷}:={\widehat{𝑷}}_j$	Estimates of $j$-th subspace
${\widehat{𝑷}}_{\mathrm{rot},k}:={\widehat{𝑷}}_{j,\mathrm{rot},k}$	$k$-th estimate of ${𝑷}_{\mathrm{rot}}$.
${𝒆}_t:={𝑰}_{{𝒯}_t}{\left({𝚿}_{{𝒯}_t}{}_{}{}^{\prime }𝚿_{{𝒯}_t}^{}\right)}^{-1}{𝑰}_{{𝒯}_t}{}_{}{}^{\prime }𝚿\left({\ell }_t+{𝒗}_t\right)$ (Proved in Lemma 5.15)	Expression for error, ${𝒆}_t={\widehat{𝒙}}_t-{𝒙}_t={\ell }_t-\hat{\ell }+{𝒗}_t$
Scalars (Derived in Fact 5.14)
${\zeta }_0^{+}:=\tilde{\epsilon }+\left|\sin \theta \right|$	Bound on $\mathrm{SE}(\widehat{𝑷},{𝑷}_{\ast })$, i.e., before subspace update
${\zeta }_1^{+}:=0.4\cdot 1.2\left(\left(0.1+\tilde{\epsilon }\right)\left|\sin \theta \right|+\tilde{\epsilon }\right)+0.11\tilde{\epsilon }$	Bound on $\mathrm{SE}([{\widehat{𝑷}}_{\ast },{\widehat{𝑷}}_{\mathrm{rot},1}],{𝑷}_{\mathrm{rot}})$, i.e., after $1$st subspace update
${\zeta }_k^{+}:=0.4\cdot \left(1.2{\zeta }_{k-1}^{+}\right)+0.11\tilde{\epsilon }$	Bound on $\mathrm{SE}([{\widehat{𝑷}}_{\ast },{\widehat{𝑷}}_{\mathrm{rot},k}],{𝑷}_{\mathrm{rot}})$, i.e., after $k$-th subspace update
Events – $t_j$ known (Proved in Sec. V-D)
${\Gamma }_0:=\left\{\mathrm{SE}({\widehat{𝑷}}_{\ast },{𝑷}_{\ast })\le \tilde{\epsilon }\right\}$	Previous subspace, ${𝑷}_{\ast }$ is $\tilde{\epsilon }$-accurately estimated
${\Gamma }_k:={\Gamma }_{k-1}\cap \left\{\mathrm{SE}([{\widehat{𝑷}}_{\ast },{\widehat{𝑷}}_{\mathrm{rot},k}],{𝑷}_{\mathrm{rot}})\le {\zeta }_k^{+}\right\}$	All $k$ subspace update steps work
${\Gamma }_{K+1}:={\Gamma }_K\cap \left\{\mathrm{SE}(\widehat{𝑷},𝑷)\le \tilde{\epsilon }\right\}$	Current subspace, $𝑷$ is $\tilde{\epsilon }$-accurately estimated
Events – $t_j$ unknown (This section is only used in Appendix A)
${\hat{t}}_{j-1,fin}:={\hat{t}}_{j-1}+\left(K+1\right)\alpha -1$	Time at which deletion step is complete
$t_{j,\ast }={\hat{t}}_{j-1,fin}+\lceil \frac{t_j-{\hat{t}}_{j-1,fin}}{\alpha }\rceil \alpha$	First possible time instant at which subspace change can be detected
$\mathrm{Det0}:=\left\{{\hat{t}}_j=t_{j,\ast }\right\}$ (Proved in Appendix A)	Subspace Change detected within $\alpha$-frames
$\mathrm{Det1}:=\left\{{\hat{t}}_j=t_{j,\ast }+\alpha \right\}$ (Proved in Appendix A)	Subspace Change detected after $\alpha$, but before $2\alpha$ frames
${\mathrm{ProjSVD}}_k:=\left\{\mathrm{SE}\left([{\widehat{𝑷}}_{\ast },{\widehat{𝑷}}_{\mathrm{rot},k}]\right)\le {\zeta }_{\mathrm{rot},k}^{+}\right\}$ (Proved in Lemma 5.16)	$k$-th Proj SVD works
$\mathrm{ProjSVD}:={\cap }_{k=1}^K{\mathrm{ProjSVD}}_k$	All $K$ Proj SVD steps work
$\mathrm{Del}:=\left\{\mathrm{SE}(\widehat{𝑷},𝑷)\le \tilde{\epsilon }\right\}$ (Proved in Lemma 5.17)	Deletion Step works
$\mathrm{NoFalseDets}$ (Proved in Appendix A)	No false detection of subspace change
${\Gamma }_{0,\mathrm{end}}:=\left\{\mathrm{SE}({\widehat{𝑷}}_{\ast },{𝑷}_{\ast })\le \tilde{\epsilon }\right\}$ (Proved in Claim 5.13)	Previous subspace estimated within $\tilde{\epsilon }$-accuracy
${\Gamma }_{j,\mathrm{end}}$ (Proved in Claim 5.13)	All previous $j$ subspaces estimated within $\tilde{\epsilon }$-accuracy
Notation for PCA in data-dependent noise: Theorem 5.7 (Proved in Appendix B)
${𝒚}_t={\ell }_t+{𝒘}_t+{𝒛}_t$	Observations: True data - ${\ell }_t=𝑷{𝒂}_t=[{𝑷}_{\mathrm{fix}},{𝑷}_{\mathrm{rot}}]{𝒂}_t$;
	Data-dep noise - ${𝒘}_t$; Modeling error - ${𝒛}_t$
${𝒘}_t={𝑴}_t{\ell }_t={𝑴}_{2,t}{𝑴}_{1,t}{\ell }_t$	Data-dependent noise,
$\Vert {𝑴}_{1,t}{𝑷}_{\ast }\Vert \le q_0$ and $\Vert {𝑴}_{1,t}{𝑷}_{\mathrm{rot}}\Vert \le q_{\mathrm{rot}}$	Assumptions on Data-dependency matrices
$\Vert {𝑴}_{2,t}\Vert \le 1$ and $\parallel \frac{1}{\alpha }{\sum }_t{𝑴}_{2,t}{𝑴}_{2,t}{}_{}{}^{\prime }\parallel \le b_0=0.01$
$\Vert {𝒛}_t\Vert \le b_z:=q_0\sqrt{r{\lambda }^{+}}+q_{\mathrm{rot}}\sqrt{{\lambda }_{\mathrm{ch}}}$, $\parallel 𝔼[{𝒛}_t{𝒛}_t\parallel \le {\lambda }_z^{+}:=b_z^2/r$	Assumptions on modeling error.

Table 5. TABLE V: Average subspace error SE ​ ( 𝑷 ^ ( t ) , 𝑷 ( t ) ) SE subscript ^ 𝑷 𝑡 subscript 𝑷 𝑡 \mathrm{SE}(\hat{\bm{P}}_{(t)},\ \bm{P}_{(t)}) and time comparison for different values of signal size n 𝑛 n . The values in brackets denote average time taken per frame (– indicates that the algorithm does not work).

	ReProCS	GRASTA	ORPCA	Offline ReProCS	PCP	AltProj	RPCA-GD
$n=500$ (in ${10}^{-4}s$)	$0.066$ ($3.1$)	$0.996$ ($2.8$)	$0.320$ ($10$)	$8.25\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}$ ($6.3$)	$1.00$ ($51$)	$0.176$ ($104$)	$0.215$ ($454$)
$n=500$, Bern. (in ${10}^{-4}s$)	$0.044$ ($4.8$)	$0.747$ ($1.9$)	$0.078$ ($1.8$)	$3.9\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟕}}$ ($9.2$)	$1.2\times {10}^{-4}$ ($395$)	$0.0001$ ($32$)	$0.303$ ($329$)
$n=5000$ (in ${10}^{-2}s$)	$0.048$ ($3.7$)	$0.999$ ($0.11$)	$0.322$ ($0.30$)	$6.05\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}$ ($8.5$)	$0.999$ ($8.9$)	$0.354$ ($13.0$)	$0.223$ ($47.0$)
$n=10,000$ (in ${10}^{-2}s$)	$0.090$ ($15.6$)	$0.999$ ($0.25$)	$0.3235$ ($0.68$)	$0.0006$ ($36.8$)	–	–	–

Equations

$$\mathrm{SE}(\hat{\mathbf{P}},\mathbf{P}):=\Vert (\mathbf{I}-\hat{\mathbf{P}}{\hat{\mathbf{P}}}^{\prime })\mathbf{P}\Vert$$

$${\mathbf{y}}_t:={\mathbf{\ell }}_t+{\mathbf{x}}_t+{\mathbf{v}}_t,\text{ for }t=1,2,\dots ,d$$

$${\mathbf{P}}_{(t)}={\mathbf{P}}_{(t_j)}\text{ for all }t\in [t_j,t_{j+1}),j=0,1,\dots ,J,$$

$$\mathrm{SE}({\mathbf{P}}_{j-1},{\mathbf{P}}_j)=\mathrm{SE}({\mathbf{P}}_{j-1,\mathrm{ch}},{\mathbf{P}}_{j,\mathrm{rot}})$$

$$\mathrm{SE}({\mathbf{P}}_{j-1},{\mathbf{P}}_j)=\mathrm{SE}({\mathbf{P}}_{j-1,\mathrm{ch}},{\mathbf{P}}_{j,\mathrm{rot}})\le \Delta .$$

$${\mathbf{P}}_{j,\mathrm{new}}:=\frac{(\mathbf{I}-{\mathbf{P}}_{j-1,\mathrm{ch}}{\mathbf{P}}_{j-1,\mathrm{ch}}{}^{\prime }){\mathbf{P}}_{j,\mathrm{rot}}}{\mathrm{SE}({\mathbf{P}}_{j-1,\mathrm{ch}},{\mathbf{P}}_{j,\mathrm{rot}})}$$

$${\mathbf{P}}_j=[\underset{{\mathbf{P}}_{j-1,\mathrm{fix}}}{\underbrace{({\mathbf{P}}_{j-1}{\mathbf{U}}_j{)}_{[1,r-1]}}},{\mathbf{P}}_{j,\mathrm{rot}}],\text{ where }$$

$${\mathbf{P}}_{j,\mathrm{rot}}:=\underset{{\mathbf{P}}_{j-1,\mathrm{ch}}}{\underbrace{({\mathbf{P}}_{j-1}{\mathbf{U}}_j{)}_r}}\cos {\theta }_j+{\mathbf{P}}_{j,\mathrm{new}}\sin {\theta }_j$$

$$\underset{j=1,2,\dots ,J}{\max }\underset{i=1,2,\dots ,n}{\max }\Vert \mathrm{basis}([{\mathbf{P}}_{j-1},{\mathbf{P}}_j]{)}^i\Vert \le \sqrt{\frac{\mu (r+1)}{n}}$$

$$\underset{j=1,2,\dots r}{\max }\underset{t}{\max }\frac{({\mathbf{a}}_t{)}_j^2}{{\lambda }_j(\mathbf{\Lambda })}\le \eta .$$

$$\gamma (\mathcal{J}):=\underset{i=1,2,\dots ,n}{\max }\frac{1}{|\mathcal{J}|}\sum _{t\in \mathcal{J}}\text{mathbbm}{1}_{\{i\in {\mathcal{T}}_t\}}.$$

$${\text{max-outlier-frac-row}}^{\alpha }:=\underset{{\mathcal{J}}^{\alpha }\subseteq [t_1,d]}{\max }\gamma ({\mathcal{J}}^{\alpha }).$$

$$K:=\lceil C\log (\Delta /\tilde{\epsilon })\rceil ,\text{ and }\alpha \ge {\alpha }_{\ast }:=C{f}^2(r\log n).$$

$$\left\{\begin{array}[]{ll}2\tilde{\varepsilon}+\Delta&\text{ if }t\in[t_{j},{\hat{t}}_{j}+\alpha)\\ 1.2\tilde{\varepsilon}+(0.5)^{k-2}0.06\Delta&\text{ if }t\in[{\hat{t}}_{j}+(k-1)\alpha,{\hat{t}}_{j}+k\alpha)\\ 2\tilde{\varepsilon}&\text{ if }t\in[{\hat{t}}_{j}+K\alpha,{\hat{t}}_{j}+K\alpha+\alpha)\\ \tilde{\varepsilon}&\text{ if }t\in[{\hat{t}}_{j}+K\alpha+\alpha,t_{j+1})\end{array}\right.$$

$$\mathrm{SE}({\hat{\mathbf{P}}}_{(t)}^{\mathrm{offline}},{\mathbf{P}}_{(t)})\le 2\tilde{\epsilon },\text{ and }\Vert {\hat{\mathbf{\ell }}}_t^{\mathrm{offline}}-{\mathbf{\ell }}_t\Vert \le 2.4\tilde{\epsilon }\Vert {\mathbf{\ell }}_t\Vert .$$

$$c_1\Delta (c_{2}f\sqrt{{\rho }_{\mathrm{row}}}{)}^K+0.2\tilde{\epsilon }\le \tilde{\epsilon }.$$

$${\mathbf{e}}_t$$

$$\mathbf{P}=[{\mathbf{P}}_{\mathrm{fix}},{\mathbf{P}}_{\mathrm{rot}}],\text{ where }{\mathbf{P}}_{\mathrm{rot}}:=({\mathbf{P}}_{\mathrm{ch}}\cos \theta +{\mathbf{P}}_{\mathrm{new}}\sin \theta ),$$

$${\mathbf{D}}_{obs}:=\frac{1}{\alpha }\sum _{t\in {\mathcal{J}}^{\alpha }}\mathbf{\Phi }{\mathbf{y}}_t{\mathbf{y}}_t{}^{\prime }\mathbf{\Phi }.$$

$${\alpha }_0:=C\eta \max (f(r\log n),\eta f^2(r+\log n)).$$

$$\Vert \frac{1}{\alpha }\sum _{t\in {\mathcal{J}}^{\alpha }}{\mathbf{M}}_{2,t}{\mathbf{M}}_{2,t}{}^{\prime }\Vert \le b_0=\mathrm{0.01.}$$

$$q_0\le 2\tilde{\epsilon },q_{\mathrm{rot}}\le 0.2|\sin \theta |,\tilde{\epsilon }f\le 0.01|\sin \theta |,\text{ and }$$

$$b_z\le C(q_0\sqrt{r{\lambda }^{+}}+q_{\mathrm{rot}}\sqrt{{\lambda }_{\mathrm{ch}}})$$

$$\mathrm{SE}(\hat{\mathbf{P}},{\mathbf{P}}_{\mathrm{rot}})$$

$$\le 1.01|\sin \theta |\frac{0.39q_{\mathrm{rot}}+0.1\tilde{\epsilon }}{|\sin \theta |}$$

$$<0.4q_{\mathrm{rot}}+0.11\tilde{\epsilon }.$$

$${\lambda }_{\max }({\mathbf{D}}_{obs})$$

$$\ge (0.97{\sin }^2\theta -0.4q_{\mathrm{rot}}|\sin \theta |-0.15\tilde{\epsilon }|\sin \theta |){\lambda }_{\mathrm{ch}}.$$

$$\Vert \frac{1}{\alpha }\sum _{t\in {\mathcal{J}}^{\alpha }}\mathbb{E}[{\mathbf{w}}_t{\mathbf{w}}_t{}^{\prime }]\Vert \le \sqrt{b_0}{c}_w\text{ and }$$

$$\Vert \frac{1}{\alpha }\sum _{t\in {\mathcal{J}}^{\alpha }}\mathbb{E}[\mathbf{\Phi }{\mathbf{\ell }}_t{\mathbf{w}}_t{}^{\prime }\mathbf{\Phi }]\Vert \le \sqrt{b_0}{c}_{wl}.$$

$$\Vert \frac{1}{\alpha }\sum _{t\in {\mathcal{J}}^{\alpha }}{\mathbf{I}}_{{\mathcal{T}}_t}{\mathbf{I}}_{{\mathcal{T}}_t}{}^{\prime }\Vert$$

$$\le {\rho }_{\mathrm{row}}=0.01/f^2\le 0.01=b_0.$$

$${\mathbf{C}}_t=\mathrm{diag}\left((c_t{)}_1,(c_t{)}_2,\cdots ,(c_t{)}_n\right),$$

Full text

I Introduction

Principal Components Analysis (PCA) is a widely used dimension reduction technique in a variety of scientific applications. Given a set of data vectors, PCA tries to finds a smaller dimensional subspace that best approximates a given dataset. According to its modern definition [2], robust PCA (RPCA) is the problem of decomposing a given data matrix into the sum of a low-rank matrix (true data) and a sparse matrix (outliers). The column space of the low-rank matrix then gives the desired principal subspace (PCA solution). In recent years, the RPCA problem has been extensively studied, e.g., [2, 3, 4, 5, 6, 7, 8]. A common application of RPCA is in video analytics in separating video into a slow-changing background image sequence (modeled as a low-rank matrix) and a foreground image sequence consisting of moving objects or people (sparse) [2]. Dynamic RPCA refers to the dynamic (time-varying) extension of RPCA [5, 9, 7]. It assumes that the true (uncorrupted) data lies in a low-dimensional subspace that can change with time, albeit slowly. This is a more appropriate model for long data sequences, e.g., surveillance videos. The goal is to track this changing subspace over time in the presence of sparse outliers. Hence this problem can also be referred to as robust subspace tracking.

I-A Notation and Problem Setting

I-A1 Notation

We use bold lower case letters to denote vectors, bold upper case letters to denote matrices, and calligraphic letters to denote sets or events. We use the interval notation $[a,b]$ to mean all of the integers between $a$ and $b$, inclusive, and $[a,b):=[a,b-1]$. We will often use $\mathcal{J}$ to denote a time interval and $\mathcal{J}^{\alpha}$ to denote a time interval of length $\alpha$. We use $\mathbbm{1}_{S}$ to denote the indicator function for statement $S$, i.e. $\mathbbm{1}_{S}=1$ if $S$ holds and $\mathbbm{1}_{S}=0$ otherwise. We use $\|\cdot\|$ without a subscript to denote the $l_{2}$ norm of a vector or the induced $l_{2}$ norm of a matrix. For other $l_{p}$ norms, we use $\|\cdot\|_{p}$. For a set $\mathcal{T}$, we use $\bm{I}_{\mathcal{T}}$ to refer to an $n\times|\mathcal{T}|$ matrix of columns of the identity matrix indexed by entries in $\mathcal{T}$. For a matrix $\bm{A}$, $\bm{A}^{\prime}$ denotes its transpose and $\bm{A}_{\mathcal{T}}:=\bm{AI}_{\mathcal{T}}$ is the sub-matrix of $\bm{A}$ that contains the columns of $\bm{A}$ indexed by entries in $\mathcal{T}$. Also, we use $\bm{A}^{i}$ to denote its $i$-th row. We use $\lambda_{\min}(.)$ ($\sigma_{\min}(.)$) to denote the minimum eigen (singular) value of a matrix. Similarly for $\lambda_{\max}(.)$ and $\sigma_{\max}(.)$. We use $\delta_{s}(\bm{A})$ to denote the $s$-restricted isometry constant (RIC) [10] of $\bm{A}$.

A matrix with mutually orthonormal columns is referred to as a basis matrix and is used to represent the subspace spanned by its columns. For basis matrices $\hat{\bm{P}}$, $\bm{P}$, we use

[TABLE]

to quantify the subspace error (SE) between their respective column spans. This measures the sine of the maximum principal angle between the subspaces. When $\hat{\bm{P}}$ and $\bm{P}$ are of the same size, then $\mathrm{SE}(.)$ is symmetric, i.e., $\mathrm{SE}(\hat{\bm{P}},\bm{P})=\mathrm{SE}(\bm{P},\hat{\bm{P}})$. We use $\bm{P}_{\perp}$ to denote a basis matrix for the orthogonal complement of $\operatorname{span}(\bm{P})$.

For a matrix $\bm{M}$, we use $\operatorname{basis}(\bm{M})$ to denote a basis matrix whose columns span the same subspace as the columns of $\bm{M}$.

The letters $c$ and $C$ denote different numerical constants in each use; $c$ is used for constants less than one and $C$ for those equal to or greater than one.

I-A2 Dynamic RPCA or Robust Subspace Tracking Problem Statement

At each time $t$, we observe $\bm{y}_{t}\in\mathbb{R}^{n}$ that satisfies

[TABLE]

where $\bm{x}_{t}$ is the sparse outlier vector, $\bm{\ell}_{t}$ is the true data vector that lies in a fixed or slowly changing low-dimensional subspace of $\mathbb{R}^{n}$, and $\bm{v}_{t}$ is small unstructured noise or modeling error. To be precise, $\bm{\ell}_{t}=\bm{P}_{(t)}\bm{a}_{t}$ where $\bm{P}_{(t)}$ is an $n\times r$ basis matrix with $r\ll n$ and with $\|(\bm{I}-\bm{P}_{(t-1)}\bm{P}_{(t-1)}{}^{\prime})\bm{P}_{(t)}\|$ small compared to $\|\bm{P}_{(t)}\|=1$ (slow subspace change). We use $\mathcal{T}_{t}$ to denote the support set of $\bm{x}_{t}$ and we let $s:=\max_{t}|\mathcal{T}_{t}|$. Given an initial subspace estimate, $\hat{\bm{P}}_{0}$, the goal is to track $\operatorname{span}(\bm{P}_{(t)})$ within a short delay of each subspace change. The initial estimate can be obtained by applying any static (batch) RPCA technique, e.g., PCP [2] or AltProj [6], to the first $t_{\mathrm{train}}$ data frames, $\bm{Y}_{[1,t_{\mathrm{train}}]}$. A by-product of our solution approach is that the true data vectors $\bm{\ell}_{t}$, the sparse outliers $\bm{x}_{t}$, and their support sets $\mathcal{T}_{t}$ can also be tracked on-the-fly. In many practical applications, in fact, $\bm{x}_{t}$ or $\mathcal{T}_{t}$ is often the quantity of interest.

We also assume that (i) $|\mathcal{T}_{t}|/n$ is upper bounded, (ii) $\mathcal{T}_{t}$ changes enough over time so that any one index is not part of the outlier support for too long, (iii) the columns of $\bm{P}_{(t)}$ are dense (non-sparse), and (iv) the subspace coefficients $\bm{a}_{t}$ are element-wise bounded, mutually independent, zero mean, have identical and diagonal covariance matrices, and are independent of the outlier supports $\mathcal{T}_{t}$. We quantify everything in Sec. II.

I-A3 Subspace Change Assumption

To ensure that the number of unknowns is not too many (see the discussion in Sec. I-C), we will further assume that the subspace $\operatorname{span}(\bm{P}_{(t)})$ is piecewise constant with time, i.e.,

[TABLE]

with $t_{0}=1$ and $t_{J+1}=d$. Let $\bm{P}_{j}:=\bm{P}_{(t_{j})}$. At each change time, $t_{j}$, the change is “slow”. This means two things:

First, at each $t_{j}$, only one direction can change with the rest of the subspace remaining fixed, i.e.,

[TABLE]

where $\bm{P}_{j-1,{\mathrm{ch}}}$ is a direction from $\operatorname{span}(\bm{P}_{j-1})$ that “changes” at $t_{j}$ and $\bm{P}_{j,{\mathrm{rot}}}$ is its “rotated” version. Thus $\operatorname{span}(\bm{P}_{j-1})=\operatorname{span}([\bm{P}_{j-1,{\mathrm{fix}}},\bm{P}_{j-1,{\mathrm{ch}}}])$ and $\operatorname{span}(\bm{P}_{j})=\operatorname{span}([\bm{P}_{j-1,{\mathrm{fix}}},\bm{P}_{j,{\mathrm{rot}}}])$ where $\bm{P}_{j-1,{\mathrm{fix}}}$ is an $n\times(r-1)$ matrix that denotes the part of the subspace that remains “fixed” at $t_{j}$.

Of course at different $t_{j}$’s, the changing directions could be different. 2. 2.

Second, the angle of change is small, i.e., for a $\Delta\ll 1$,

[TABLE]

I-A4 Equivalent generative model

With the above model,

[TABLE]

is the newly added direction at $t_{j}$, $\theta_{j}:=\cos^{-1}|\bm{P}_{j-1,{\mathrm{ch}}}{}^{\prime}\bm{P}_{j,{\mathrm{rot}}}|$ is the angle by which $\bm{P}_{j-1,{\mathrm{ch}}}$ gets rotated out-of-plane (towards $\bm{P}_{j,\mathrm{new}}$ which lies in $\operatorname{span}(\bm{P}_{j-1})^{\perp}$) to get $\bm{P}_{j,{\mathrm{rot}}}$. Without loss of generality, assume $0\leq\theta_{j}\leq\pi/2$. Thus,

• •

$|\sin\theta_{j}|=\sin\theta_{j}=\mathrm{SE}(\bm{P}_{j-1,{\mathrm{ch}}},\bm{P}_{j,{\mathrm{rot}}})=\mathrm{SE}(\bm{P}_{j-1},\bm{P}_{j})\leq\Delta$, and

• •

$\bm{P}_{j,\mathrm{del}}:=\bm{P}_{j-1,{\mathrm{ch}}}\sin\theta_{j}-\bm{P}_{j,\mathrm{new}}\cos\theta_{j}$ is the direction that got deleted at $t_{j}$.

We have the following equivalent generative model for getting $\bm{P}_{j}$ from $\bm{P}_{j-1}$: let $\bm{U}_{j}$ be an $r\times r$ rotation matrix,

[TABLE]

For a simple example of this in 3D ($n=3$), see Fig. 1.

To make our notation easy to remember, we try to explain its meaning better. Consider the change at $t_{j}$. The direction from $\operatorname{span}(\bm{P}_{j-1})$ that changes is denoted by $\bm{P}_{j-1,{\mathrm{ch}}}$. This changes by getting rotated (out-of-plane) by a small angle $\theta_{j}$ towards a new out-of-plane direction $\bm{P}_{j,\mathrm{new}}$ to get the changed/rotated direction $\bm{P}_{j,{\mathrm{rot}}}$. Here “plane” refers to the hyperplane $\operatorname{span}(\bm{P}_{j-1})$. The basis for the $r-1$-dimensional subspace of $\operatorname{span}(\bm{P}_{j-1})$ that does not change at $t_{j}$ is $\bm{P}_{j-1,{\mathrm{fix}}}$. So $\bm{P}_{j}=[\bm{P}_{j-1,{\mathrm{fix}}},\bm{P}_{j,{\mathrm{rot}}}]$.

The span of left singular vectors of $\bm{L}$ is contained in, or equal to, $\operatorname{span}([\bm{P}_{0},\bm{P}_{1,\mathrm{new}},\bm{P}_{2,\mathrm{new}},\dots,\bm{P}_{J,\mathrm{new}}])$. Equality holds if $\bm{P}_{j,\mathrm{new}}$ is orthogonal to $\operatorname{span}([\bm{P}_{0},\bm{P}_{1,\mathrm{new}},\dots,\bm{P}_{j-1,\mathrm{new}}])$ for each $j$.

In this work we have assumed the simplest possible model on subspace change where, at a change time, only one direction can change. Observe though that, at different change times, the changing direction could be different and hence, over a long period of time, the entire subspace could change. This simple model can be generalized to $r_{\mathrm{ch}}>1$ directions changing; see the last appendix in the ArXiv posting of this work. It is also possible to study the most general case where $r_{\mathrm{ch}}=r$ and hence no model is assumed for subspace change (only a bound on the maximum principal angle of the change). This requires significant changes to both the algorithm and the guarantee; it is studied in follow-up work [11].

I-A5 Relation to original RPCA

To connect with the original RPCA problem [2, 4, 6], define the $n\times d$ data matrix $\bm{Y}:=[\bm{y}_{1},\bm{y}_{2},\dots\bm{y}_{d}]:=\bm{L}+\bm{X}+\bm{V}$ where $\bm{L}$, $\bm{X}$, $\bm{V}$ are similarly defined. Let $r_{\scriptscriptstyle{L}}$ denote the rank of $\bm{L}$ and use max-outlier-frac-col and max-outlier-frac-row to denote the maximum fraction of outliers per column and per row of $\bm{Y}$. RPCA results bound $\max(\text{\small{max-outlier-frac-row}},\text{\small{max-outlier-frac-col}})$. For dynamic RPCA, we will define max-outlier-frac-row slightly differently. It will be the maximum fraction per row of any $n\times\alpha$ sub-matrix of $\bm{Y}$ with $\alpha$ consecutive columns. Here $\alpha$ denotes the number of frames used in each subspace update. We will denote this by $\text{\small{max-outlier-frac-row}}^{\alpha}$ to indicate the difference. Since $\alpha$ is large enough (see (9)), the two definitions are only a little different. The dynamic RPCA assumption of a bound on $\max_{t}|\mathcal{T}_{t}|/n$ is equivalent to bounding max-outlier-frac-col since $\text{\small{max-outlier-frac-col}}=\max_{t}|\mathcal{T}_{t}|/n$. The requirement of $\mathcal{T}_{t}$’s changing enough is equivalent to a bound on $\text{\small{max-outlier-frac-row}}^{\alpha}$. As we explain later, the denseness assumption on the $\bm{P}_{(t)}$’s is similar to the denseness (incoherence) of left singular vectors of $\bm{L}$ assumed by all standard RPCA solutions, while the assumptions on $\bm{a}_{t}$’s replace the right singular vectors’ incoherence assumption of standard RPCA.

I-B Related Work and our Contributions

I-B1 Related Work

We briefly mention all related work here, but provide a detailed discussion later in Sec. III. There is very little work on other solutions for provably correct dynamic RPCA. This includes our early work on a partial guarantee (guarantee required assumptions on intermediate algorithm estimates) [5] and later complete correctness results [12, 7] for more complicated ReProCS-based algorithms. We refer to all of these as “original-ReProCS”. It also includes older work on modified-PCP, which is a batch solution for RPCA with partial subspace knowledge, and which can be shown to also provably solve dynamic RPCA in a piecewise batch fashion [13]. The original-ReProCS guarantees require strong assumptions on how the outlier support changes (need a very specific model inspired by a video moving object); their subspace change assumptions are unrealistic; and their subspace tracking delay (equal to the required delay between subspace change times) is very large. On the other hand, the Modified-PCP guarantee [13] requires the outlier support to be uniformly randomly generated (strong assumption; for video, it means that the moving objects need to be single pixel wide and should be jumping around randomly from frame to frame); requires a different stronger assumption on subspace change; and cannot detect subspace change automatically. Other than the above, there is some work on online algorithms for RPCA. The only work that comes with some guarantee, although it is a partial guarantee, is an online stochastic optimization based solver for the PCP convex program (ORPCA) [14]. Its guarantee assumed that the basis matrix for the subspace estimate at each $t$ was full rank. To our best knowledge, there is no follow-up work on a complete correctness result for it. There is also much work on empirical online solutions for RPCA, e.g., [15], and older work, e.g., [16, 17]. From a practical standpoint, any online algorithm will implicitly also provide a tracking solution. However, as shown in Sec. VI, the solution is not as good as that of ReProCS which explicitly exploits slow subspace change.

The standard RPCA problem has been extensively studied [2, 3, 4, 6, 8, 18]. We discuss these works in detail in Sec. III. A summary is provided in Table I. Briefly, these either need outlier fractions in each row and each column of the observed data matrix to be $O(1/r_{\scriptscriptstyle{L}})$ (AltProj [6], GD [8], NO-RMC [18], PCP result of [3, 4], denoted PCP(H)) or need the outlier support to be uniformly randomly generated (PCP result of [2], denoted PCP(C)). Moreover, all these are batch solutions with large memory complexity $O(nd)$.

I-B2 Contributions

We develop a simple algorithm, termed simple-ReProCS or s-ReProCS, for provably solving the robust subspace tracking or dynamic RPCA problem described earlier. We also develop its offline extension that can be directly compared with the standard RPCA results. Simple-ReProCS is based on the ReProCS framework [5]. Our main contribution is the first correctness guarantee for dynamic RPCA that holds under weakened versions of standard RPCA assumptions, slow subspace change, and a lower bound on most outlier magnitudes (this lower bound is proportional to the rate of subspace change). We say “weakened” because our guarantee implies that, after initialization, s-ReProCS can tolerate an order-wise larger fraction of outliers per row than all existing approaches, without requiring the outlier support to be uniformly randomly generated or without needing any other model on support change. It allows $\text{\small{max-outlier-frac-row}}^{\alpha}\in O(1)$ (instead of $O(1/r_{\scriptscriptstyle{L}})$). For the video application, this implies that it tolerates slow moving and occasionally static foreground objects much better than other approaches. This fact is also backed up by comparisons on real videos, see Sec. VI and also see [21].

A second key contribution is the algorithm itself. Unlike original-ReProCS [12, 7], s-ReProCS ensures that the estimated subspace dimension is bounded by $(r+1)$ at all times without needing the complicated cluster-EVD step. More importantly, s-ReProCS is provably fast and memory-efficient: its time complexity is comparable to that of SVD for vanilla PCA, and its memory complexity is near-optimal and equal to $O(nr\log n\log(1/\tilde{\varepsilon}))$ where $\tilde{\varepsilon}$ is the desired subspace recovery accuracy. This is near-optimal because $nr$ is the memory needed to output an $r$-dimensional subspace estimate in $\mathbb{R}^{n}$, and the complexity is within log factors of the optimal. To our best knowledge, s-ReProCS is the first provably correct RPCA or dynamic RPCA solution that is as fast as the best RPCA solution in terms of computational complexity without requiring the data matrix to be nearly square and has near-optimal memory complexity. We provide a tabular comparison of guarantees of offline s-ReProCS with other provable RPCA solutions in Table I. We compare s-ReProCS with other online or tracking solutions for RPCA or dynamic RPCA in Table II (original-ReProCS, modified-PCP, follow-up work on RePrOCS-NORST [19, 11], ORPCA and GRASTA).

We give a significantly shorter and simpler proof than that for the earlier guarantees for ReProCS-based methods. We do this by first separately proving a result for the problem of “correlated-PCA” or “PCA in data-dependent noise” [22, 23] with partial subspace knowledge. This result given in Theorem 5.7 of Sec. V-A may also be of independent interest.

I-C The need for a piecewise constant model on subspace change

We explain why the piecewise-constant subspace change model is needed. Even if the observed data were perfect (no noise/outlier/missing-data, i.e., we observed $\bm{\ell}_{t}$, and all measurements were linearly independent) and the previous subspace were exactly known, in order to obtain a correct $r$-dimensional estimate111requires finding both the newly added direction, $\bm{P}_{j,\mathrm{new}}$, and the deleted direction, $\bm{P}_{j,\mathrm{del}}$ for each $\bm{P}_{(t)}$, one would need at least $r$ samples. Of course, to just find the newly added direction $\bm{P}_{j,\mathrm{new}}$ and use an $(r+1)$-dimensional estimate, one sample would suffice in this ideal setting (doing this will be especially problematic if the subspace changes at each time because it will mean the estimated subspace dimension will keep growing as $r+t$ at time $t$). Our actual setting is not this ideal one: we know the previous subspace only up to $\epsilon$ error and we observe $\bm{y}_{t}$ which is a noisy and outlier-corrupted version of $\bm{\ell}_{t}$. This is why, in our setting, more than one data samples are needed even to accurately estimate the newly added direction. Since we get only one observed data vector $\bm{y}_{t}$ at each time, the only way to have enough data samples for estimating each subspace is to assume that $\bm{P}_{(t)}$ is piecewise constant with time, i.e., it satisfies (2). In fact, our required lower bound on $t_{j+1}-t_{j}$ is only a little more than $r$ (see Theorem 2.2), thus making our model a good approximation to slow continuous subspace change.

Furthermore, the following point should be mentioned. In the entire literature on subspace tracking (both with and without outliers, and with and without even missing data), there is no model for subspace change for which there are any provable guarantees. There is no work on provable subspace tracking with outliers (robust subspace tracking) except our own previous work which also used the piecewise constant subspace change model. The subspace tracking (ST) problem (without outliers), and with or without missing data, has been extensively studied [24, 25, 26, 27, 28, 29]; however, all existing guarantees are asymptotic results for the statistically stationary setting of data being generated from a single unknown subspace. Moreover, most of these also make assumptions on intermediate algorithm estimates. For a longer discussion of this, please see [30].

I-D Paper Organization

The proposed algorithm, simple-ReProCS, and its performance guarantees, Theorem 2.2, are given in Sec. II. We discuss the related work in detail in Sec. III and explain how our guarantee compares with other provable results on RPCA or dynamic RPCA from the literature. Sec. IV provides the main ideas that lead to the proof of Theorem 2.2. We prove Theorem 2.2 under the assumption that the subspace change times are known in Sec. V. This proof helps illustrate all the ideas of the actual proof but with minimal notation. The general proof of Theorem 2.2 is given in Appendix A. Theorem 2.2 relies on a guarantee for PCA in data-dependent noise [22, 23] when partial subspace knowledge is available. This result is proved in Appendix B. We provide detailed empirical evaluation evaluation of simple-ReProCS in Sec. VI. We conclude and discuss future directions in Sec. VII.

II The simple-ReProCS Algorithm and its Guarantee

II-A Simple-ReProCS (s-ReProCS)

S-ReProCS proceeds as follows. The initial subspace is assumed to be accurately known (obtained using AltProj or PCP). At time $t$, if the previous subspace estimate, $\hat{\bm{P}}_{(t-1)}$, is accurate enough, because of slow subspace change, projecting $\bm{y}_{t}=\bm{x}_{t}+\bm{\ell}_{t}+\bm{v}_{t}$ onto its orthogonal complement will nullify most of $\bm{\ell}_{t}$. Moreover, $\|\bm{v}_{t}\|$ is small (by assumption). We compute $\tilde{\bm{y}}_{t}:=\bm{\Psi}\bm{y}_{t}$ where $\bm{\Psi}:=\bm{I}-\hat{\bm{P}}_{(t-1)}\hat{\bm{P}}_{(t-1)}{}^{\prime}$. Thus, $\tilde{\bm{y}}_{t}=\bm{\Psi}\bm{x}_{t}+\bm{b}_{t}$ where $\bm{b}_{t}:=\bm{\Psi}(\bm{\ell}_{t}+\bm{v}_{t})$ and $\|\bm{b}_{t}\|$ is small. Recovering $\bm{x}_{t}$ from $\tilde{\bm{y}}_{t}$ is thus a traditional compressive sensing (CS) / sparse recovery problem in small noise [10]. This is solvable because incoherence (denseness) of $\bm{P}_{(t)}$’s and slow subspace change implies [5] that $\bm{\Psi}$ satisfies the restricted isometry property [10]. We compute $\bm{\hat{x}}_{t,cs}$ using $l_{1}$ minimization followed by thresholding based support estimation to get $\hat{\mathcal{T}}_{t}$. A Least Squares (LS) based debiasing step on $\hat{\mathcal{T}}_{t}$ returns the final $\bm{\hat{x}}_{t}$. We then estimate $\bm{\ell}_{t}$ as $\hat{\bm{\ell}}_{t}=\bm{y}_{t}-\hat{\bm{x}}_{t}$. We refer to the above step as Projected Compressive Sensing (CS). As explained in [21, 30], this can also be understood as solving a Robust Regression problem222The above step equivalently solves for $\tilde{\bm{a}},\tilde{\bm{x}}$ that satisfy $\bm{y}_{t}=\hat{\bm{P}}_{t-1}\tilde{\bm{a}}+\tilde{\bm{x}}+\bm{b}_{t}$ with $\tilde{\bm{x}}$ being sparse and $\|\bm{b}_{t}\|$ being small. This is the approximate robust regression problem where columns of $\hat{\bm{P}}_{t-1}$ are the regressors/predictors, $\tilde{\bm{x}}$ is the sparse outliers and $\bm{b}_{t}$ is the small “noise” or model inaccuracy..

The $\hat{\bm{\ell}}_{t}$’s are used for the Subspace Update step which involves (i) detecting subspace change; (ii) obtaining improved estimates of the changed direction(s) by $K$ steps of projection-SVD [5], each done with a new set of $\alpha$ frames of $\hat{\bm{\ell}}_{t}$; and (iii) a simple SVD based subspace re-estimation step, done with another new set of $\alpha$ frames. This is done to remove the deleted direction and get an $r$-dimensional estimate of the new subspace. We explain the subspace change detection strategy in Sec. IV-B. Suppose the change is detected at ${\hat{t}}_{j}$. The $k$-th projection-SVD step involves computing $\hat{\bm{P}}_{j,{\mathrm{rot}},k}$ as the top singular vector of $(\bm{I}-\hat{\bm{P}}_{j-1}\hat{\bm{P}}_{j-1}{}^{\prime})[\hat{\bm{\ell}}_{{\hat{t}}_{j}+(k-1)\alpha},\hat{\bm{\ell}}_{{\hat{t}}_{j}+(k-1)\alpha+1},\dots,\hat{\bm{\ell}}_{{\hat{t}}_{j}+k\alpha-1}]$ and setting $\hat{\bm{P}}_{(t)}=\hat{\bm{P}}_{j,k}:=[\hat{\bm{P}}_{j-1},\hat{\bm{P}}_{j,{\mathrm{rot}},k}]$. For ease of understanding, we summarize a basic version of s-ReProCS in Algorithm 1. This assumes that the change times $t_{j}$ are known, i.e., that ${\hat{t}}_{j}=t_{j}$. The actual algorithm that detects changes automatically is longer and is given as Algorithm 4 in Sec. IV-B. We both analyze and implement this one.

The above approach works because, every time the subspace changes, with high probability (whp), the change can be detected within a short delay, and after that, the $K$ projection-SVD steps help get progressively improved estimates of the changed/rotated direction $\bm{P}_{j,{\mathrm{rot}}}$. The final simple SVD step re-estimates the entire subspace in order to delete $\bm{P}_{j,\mathrm{del}}$, from the estimate.

The estimates of the subspace or of $\bm{\ell}_{t}$’s are improved in offline mode as follows. At $t={\hat{t}}_{j}+K\alpha$, the $K$ projection-SVD steps are complete and hence the subspace estimate at this time is accurate enough whp. At this time, offline s-ReProCS (last line of Algorithm 4) goes back and sets $\hat{\bm{P}}_{(t)}\leftarrow[\hat{\bm{P}}_{j-1},\hat{\bm{P}}_{j,{\mathrm{rot}},K}]$ for all $t\in[{\hat{t}}_{j-1}+K\alpha,{\hat{t}}_{j}+K\alpha)$. It also uses this to get improved estimates of $\bm{\hat{x}}_{t}$ and $\hat{\bm{\ell}}_{t}$ for all these times $t$.

II-B Assumptions and Main Result

II-B1 Incoherence (denseness) of columns of $\bm{P}_{j}$’s

In order to separate the $\bm{\ell}_{t}$’s from the sparse outliers $\bm{x}_{t}$, we need an assumption that ensures that the $\bm{\ell}_{t}$’s are themselves not sparse. One way to ensure this is to assume $\mu$-incoherence [2] of the basis matrix for the subspace spanned by the columns of $\bm{P}_{j-1}$ and $\bm{P}_{j}$, i.e., assume that

[TABLE]

for a $\mu\geq 1$ but not too large (assumed to be a numerical constant henceforth). Because of our subspace change model, the subspace spanned by the columns of $[\bm{P}_{j-1},\bm{P}_{j}]$ has dimension $r+1$. In fact, $\operatorname{basis}([\bm{P}_{j-1},\bm{P}_{j}])=[\bm{P}_{j-1},\bm{P}_{j,\mathrm{new}}]$.

It is easy to see that (6), along with the bound on max-outlier-frac-col assumed in Theorem 2.2 given below ($\text{\small{max-outlier-frac-col}}\leq{0.01}/{(2\mu(r+1))}$), implies that (13) given later holds333This is true because (i) for any basis matrix $\bm{P}$, $\max_{\mathcal{T}:|\mathcal{T}|\leq 2s}\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}\|^{2}\leq 2s\max_{i}\|\bm{I}_{i}{}^{\prime}\bm{P}\|^{2}$ [5], here $s=\text{\small{max-outlier-frac-col}}\cdot n$; (ii) if $\tilde{\bm{P}}$ is such that $\operatorname{span}(\tilde{\bm{P}})\subseteq\operatorname{span}(\bm{P})$, then $\|\bm{I}_{\mathcal{T}}{}^{\prime}\tilde{\bm{P}}\|^{2}\leq\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}\|^{2}$; and (iii) both $\operatorname{span}(\bm{P}_{j})$ and $\operatorname{span}(\bm{P}_{j,\mathrm{new}})$ are contained in the span of $\operatorname{basis}([\bm{P}_{j-1},\bm{P}_{j}])$. In fact $\operatorname{basis}([\bm{P}_{j-1},\bm{P}_{j}])=[\bm{P}_{j-1},\bm{P}_{j,\mathrm{new}}]$. Thus, using the max-outlier-frac-col bound, $\max_{\mathcal{T}:|\mathcal{T}|\leq 2s}\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}_{j}\|^{2}\leq 2s\mu(r+1)/n=2\text{\small{max-outlier-frac-col}}\mu(r+1)\leq 0.01$ and the same also holds for $\max_{\mathcal{T}:|\mathcal{T}|\leq 2s}\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}_{j,\mathrm{new}}\|^{2}$.. Our result actually only needs (13), but that is complicated to state and explain. Hence we use the above stronger but well-understood assumption.

II-B2 Assumption on principal subspace coefficients $\bm{a}_{t}$

We assume that the $\bm{a}_{t}$’s are zero mean, mutually independent, element-wise bounded random variables (r.v.), have identical and diagonal covariance matrix denoted $\bm{\Lambda}$, and are independent of the outlier supports $\mathcal{T}_{t}$. Here element-wise bounded means that there exists a numerical constant $\eta$, such that

[TABLE]

For most bounded distributions, $\eta$ is a little more than one, e.g., if the entries of $\bm{a}_{t}$ are zero mean uniform, then $\eta=3$. As we explain later the above assumptions of $\bm{a}_{t}$ replace the right singular vectors’ incoherence assumption used by all standard RPCA solutions.

II-B3 Outlier fractions bounded

Similar to earlier RPCA works, we also need outlier fractions to be bounded. However, we need different bounds on this fraction per column and per row. The row bound can be much larger444One practical application where this is useful is for slow moving or occasionally static video foreground moving objects. For a stylized example of this, see Model 6.19 given in Sec. VI.. Since the ReProCS subspace update step operates on mini-batches of data of size $\alpha$ (i.e. on $n\times\alpha$ sub-matrices of consecutive columns), we need to bound max-outlier-frac-row for each such sub-matrix. We denote this by $\text{\small{max-outlier-frac-row}}^{\alpha}$.

Definition 2.1**.**

For a time interval, $\mathcal{J}$, define

[TABLE]

Thus $\gamma(\mathcal{J})$ is the maximum outlier fraction in any row of the sub-matrix $\bm{Y}_{\mathcal{J}}$ of $\bm{Y}$. Let $\mathcal{J}^{\alpha}$ denote a time interval of duration $\alpha$. Define

[TABLE] 2. 2.

Define $\text{\small{max-outlier-frac-col}}:=\max_{t}|\mathcal{T}_{t}|/n$. 3. 3.

Let $x_{\min}:=\min_{t}\min_{i\in\mathcal{T}_{t}}|(\bm{x}_{t})_{i}|$ denote the minimum outlier magnitude. 4. 4.

Use $\lambda^{-}$ and $\lambda^{+}$ to denote the minimum and maximum eigenvalues of $\bm{\Lambda}$ and $f:=\frac{\lambda^{+}}{\lambda^{-}}$ its condition number. 5. 5.

Split $\bm{a}_{t}$ as $\bm{a}_{t}=\begin{bmatrix}\bm{a}_{t,{\mathrm{fix}}}{}\\ \bm{a}_{t,{\mathrm{ch}}}\end{bmatrix}$ where $\bm{a}_{t,{\mathrm{ch}}}$ is the scalar coefficient corresponding to the changed direction. Similarly split its diagonal covariance matrix as $\bm{\Lambda}=\begin{bmatrix}\bm{\Lambda}_{\mathrm{fix}}&\bm{0}\\ \bm{0}&\lambda_{\mathrm{ch}}\end{bmatrix}$. 6. 6.

Let $\tilde{\varepsilon}$ denote the bound on initial subspace error, i.e., let $\mathrm{SE}(\hat{\bm{P}}_{0},\bm{P}_{0})\leq\tilde{\varepsilon}.$ 7. 7.

*For numerical constants $C$ that are re-used to denote different numerical values, define *

[TABLE]

II-B4 Main Result

We can now state our main result. For ease of understanding, we provide a table explaining various symbols, and assumptions required for Theorem 2.2 in Table III.

Theorem 2.2**.**

Consider simple-ReProCS given in Algorithm 4. Assume that $\mathrm{SE}(\hat{\bm{P}}_{0},\bm{P}_{0})\leq\tilde{\varepsilon}$ with $\tilde{\varepsilon}f\leq 0.01\mathrm{SE}(\bm{P}_{j-1},\bm{P}_{j})$.

(statistical assumptions) assumptions on $\bm{a}_{t}$’s hold; 2. 2.

(subspace change)

1. (a)

(2), (3), and (4) hold with $t_{j+1}-t_{j}>(K+3)\alpha$ where $K$ and $\alpha$ are defined above in (9), 2. (b)

$\Delta$* satisfies $C(2\tilde{\varepsilon}\sqrt{r\lambda^{+}}+\Delta\sqrt{\lambda_{\mathrm{ch}}})<x_{\min}/15$ with $C=\sqrt{\eta}$; * 3. 3.

(outlier fractions and left incoherence)

1. (a)

(6) holds and $\text{\small{max-outlier-frac-col}}\leq\rho_{\mathrm{col}}:=\frac{0.01}{2\mu(r+1)}$, 2. (b)

$\text{\small{max-outlier-frac-row}}^{\alpha}\leq\rho_{\mathrm{row}}:=\frac{0.01}{f^{2}}$; 4. 4.

(noise $\bm{v}_{t}$) $\bm{v}_{t}$’s are zero mean, mutually independent, independent of the $\bm{x}_{t}$’s and $\bm{\ell}_{t}$’s, and satisfy $\|\bm{v}_{t}\|^{2}\leq 0.1\tilde{\varepsilon}^{2}r\lambda^{+}$ and $\|\mathbb{E}[\bm{v}_{t}\bm{v}_{t}{}^{\prime}]\|\leq 0.1\tilde{\varepsilon}^{2}\lambda^{+}$; 5. 5.

(algorithm parameters) set $K$ and $\alpha$ as in (9), $\xi=x_{\min}/15$, $\omega_{supp}=x_{\min}/2$, $\omega_{evals}=5\tilde{\varepsilon}^{2}f\lambda^{+}$;

then, with probability at least $1-12dn^{-12}$, at all times, $t$,

$\hat{\mathcal{T}}_{t}=\mathcal{T}_{t}$, 2. 2.

$t_{j}\leq{\hat{t}}_{j}\leq t_{j}+2\alpha$, 3. 3.

$\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{(t)})\leq$* *

[TABLE] 4. 4.

and $\|\bm{\hat{x}}_{t}-\bm{x}_{t}\|=\|\hat{\bm{\ell}}_{t}-\bm{\ell}_{t}\|\leq C(\tilde{\varepsilon}\sqrt{r\lambda^{+}}+\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{(t)})\sqrt{\lambda_{\mathrm{ch}}}$ with $\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{(t)})$ bounded as above.

Consider offline s-ReProCS (last line of Algorithm 4). At all $t$,

[TABLE]

The upper bound on $\bm{v}_{t}$ and the lower bound on $x_{\min}$ can be relaxed significantly to get a more complicated result which we state in the corollary below.

Corollary 2.3**.**

Let $x_{\min,t}:=\min_{i\in\mathcal{T}_{t}}|(\bm{x}_{t})_{i}|$ denote the minimum outlier magnitude at time $t$ and define the time intervals

• •

$\mathcal{J}_{0}=[t_{j},{\hat{t}}_{j})$* (interval before the change gets detected),*

• •

$\mathcal{J}_{k}:=[{\hat{t}}_{j}+(k-1)\alpha,{\hat{t}}_{j}+k\alpha)$* ($k$-th subspace update uses data from this interval) for $k=1,2,3,\dots,K$,*

• •

and $\mathcal{J}_{K+1}:=[{\hat{t}}_{j}+K\alpha,{\hat{t}}_{j}+K\alpha+\alpha)$ (final SVD-based re-estimation step uses data from this interval).

All conclusions of Theorem 2.2 hold if the following hold instead of assumptions 1b, 4, and 5 of Theorem 2.2:

$\bm{v}_{t}$’s are zero mean, mutually independent, independent of the $\bm{x}_{t}$’s and $\bm{\ell}_{t}$’s, $\|\bm{v}_{t}\|\leq b_{v,t}$, $\|\mathbb{E}[\bm{v}_{t}\bm{v}_{t}{}^{\prime}]\|\leq b_{v,t}^{2}/r$, $x_{\min,t}$ and $b_{v,t}$ satisfy the following:

for $t\in\mathcal{J}_{0}\cup\mathcal{J}_{1}$, $b_{v,t}=C(2\tilde{\varepsilon}\sqrt{r\lambda^{+}}+0.11\Delta\sqrt{\lambda_{\mathrm{ch}}})$, and $x_{\min,t}\geq 30b_{v,t}$, 2. 2.

for $t\in\mathcal{J}_{k}$, $b_{v,t}=C(2\tilde{\varepsilon}\sqrt{r\lambda^{+}}+0.5^{k-2}0.06\Delta\sqrt{\lambda_{\mathrm{ch}}})$, and $x_{\min,t}\geq 30b_{v,t}$, for $k=2,\dots,K$, 3. 3.

for $t\in\mathcal{J}_{K+1}$, $b_{v,t}=C(\tilde{\varepsilon}\sqrt{r\lambda^{+}})$ and $x_{\min,t}\geq 30b_{v,t}$

with $C=\sqrt{\eta}$; and we set $\omega_{supp,t}=x_{\min,t}/2$, and $\xi_{t}=x_{\min,t}/15$ (alternatively, one can also set $\omega_{supp,t}$ and $\xi_{t}$ to be proportional to $b_{v,t}$ which itself is proportional to the bound on $\|\bm{x}_{t}-\bm{\hat{x}}_{t}\|$ is each interval).

Proof.

We explain the ideas leading to the proof in Sec. IV. Instead of first proving Theorem 2.2 and then Corollary 2.3, we directly only prove the latter. The proof of the former is almost the same and is immediate once the latter proof can be understood. For notational simplicity, we first prove the results under the assumption ${\hat{t}}_{j}=t_{j}$ in Sec. V. The proof without assuming ${\hat{t}}_{j}=t_{j}$ is given in Appendix A.∎

With the above corollary, the following remark is immediate.

Remark 2.4** (Bi-level outliers).**

The lower bound on outlier magnitudes can be relaxed to the following which only requires that most outlier magnitudes are lower bounded, while the others have small enough magnitudes so that their squared sum is upper bounded: Assume that the outlier magnitudes are such that the following holds: $\bm{x}_{t}$ can be split as $\bm{x}_{t}=(\bm{x}_{t})_{small}+(\bm{x}_{t})_{large}$ with the two components having disjoint supports and being such that, $\|(\bm{x}_{t})_{small}\|\leq b_{v,t}$ and the smallest nonzero entry of $(\bm{x}_{t})_{large}$ is greater than $30b_{v,t}$ with $b_{v,t}$ as defined in Corollary 2.3. If the above is true, and if the vectors $(\bm{x}_{t})_{small}$ are zero mean, mutually independent, and independent of $\bm{\ell}_{t}$’s and of the support of $(\bm{x}_{t})_{large}$, then all conclusions of Theorem 2.2 hold except the exact support recovery conclusion (this gets replaced by exact recovery of the support of $(\bm{x}_{t})_{large}$).

This remark follows by replacing $\bm{v}_{t}$ by $\bm{v}_{t}+(\bm{x}_{t})_{small}$ and $\bm{x}_{t}$ by $(\bm{x}_{t})_{large}$ in Corollary 2.3.

Remark 2.5**.**

The first condition (accurate initial estimate) can be satisfied by applying any standard RPCA solution, e.g., PCP, AltProj, or GD, on the first $t_{\mathrm{train}}=Cr$ data frames. This requires assuming that $t_{1}\geq Cr$, and that $\bm{Y}_{[1,t_{\mathrm{train}}]}$ has outlier fractions in any row or column bounded by $c/r$. Moreover, it is possible to significantly relax the initial estimate requirement to only requiring that $\mathrm{SE}(\hat{\bm{P}}_{0},\bm{P}_{0})\leq c/\sqrt{r}$ if we use $K$ iterations of the approach of follow-up work [11] to improve the estimate of $\bm{P}_{0}$ until a $\tilde{\varepsilon}$ accurate estimate is obtained, and then run s-ReProCS. For this to work, we will need a larger lower bound on $x_{\min}$ for the initial period.

Remark 2.6** (Connecting to the left incoherence of standard RPCA solutions).**

With minor changes, our left incoherence assumption, (6), can be replaced by something that is very close to the one used by all standard RPCA solutions. Instead of (6), we can assume $\mu$-incoherence of $\operatorname{basis}([\bm{P}_{0},\bm{P}_{1,\mathrm{new}},\bm{P}_{2,\mathrm{new}},\dots,\bm{P}_{J,\mathrm{new}}])$. This implies that555follows because the union of the spans of $\bm{P}_{j-1}$ and $\bm{P}_{j}$ is contained in the span of $[\bm{P}_{0},\bm{P}_{1,\mathrm{new}},\bm{P}_{2,\mathrm{new}},\dots,\bm{P}_{J,\mathrm{new}}]$. the RHS of (6) is bounded by $\sqrt{\mu(r+J)/n}$. With this, the only change to Theorem 2.2 will be that we will need $\text{\small{max-outlier-frac-col}}\leq 0.01/(2\mu(r+J))$.

If $\bm{P}_{j,\mathrm{new}}$ is orthogonal to $\operatorname{span}([\bm{P}_{0},\bm{P}_{1,\mathrm{new}},\dots,\bm{P}_{j-1,\mathrm{new}}])$ for each $j$, then the matrix $[\bm{P}_{0},\bm{P}_{1,\mathrm{new}},\bm{P}_{2,\mathrm{new}},\dots,\bm{P}_{J,\mathrm{new}}]$ is itself a basis matrix, its span is equal to that of the left singular vectors of $\bm{L}$, and its rank $r+J=r_{\scriptscriptstyle{L}}$. In this case, the above assumption and the corresponding required bound on max-outlier-frac-col are exactly the same as those used by the standard RPCA solutions.

II-C Discussion

In this section, we discuss the various implications of our result, the speed and memory guarantees, explain how to set algorithm parameters, and finally discuss its limitations.

II-C1 Subspace change detection and tracking with short delay

Theorem 2.2 shows that, whp, the subspace change gets detected within a delay of at most $2\alpha=Cf^{2}(r\log n)$ frames, and the subspace gets estimated accurately within at most $(K+3)\alpha=Cf^{2}(r\log n)\log(1/\tilde{\varepsilon})$ frames. Each column of the low rank matrix is recovered with a small time-invariant bound without any delay. If offline processing is allowed, with a delay of at most $(K+3)\alpha$, we can guarantee all recoveries within normalized error $\tilde{\varepsilon}$, or, in fact, with minor modifications, within any $\epsilon=c\tilde{\varepsilon}$ for $c<1$ (also see the limitations’ discussion). Notice that the required delay between subspace change times is more than $r$ by only logarithmic factors (assuming $f$ does not grow with $n$ or $r$). Since the previous subspace is not exactly known (is known within error at most $\tilde{\varepsilon}$), at each update step, we do need to estimate an $r$-dimensional subspace, and not a one-dimensional one. Hence it is not clear if the required delay can be reduced any further. Moreover, the delay required for the deletion step cannot be less than $r$ even in the ideal case when $\bm{\ell}_{t}$ is directly observed.

II-C2 Bi-level outliers

Consider the upper bound on $\Delta$ (amount of subspace change). Observe that the upper bound essentially depends on the ratio between $x_{\min}$ (minimum outlier magnitude) and $\sqrt{\lambda_{\mathrm{ch}}}$. Read another way, this means that $x_{\min}$ needs to be lower bounded. On first glance, this may seem counter-intuitive since sufficiently small magnitude corruptions should not be problematic. This is actually true. Sufficiently small magnitude corruptions get classified as the small noise $\bm{v}_{t}$. Moreover, as noted in Corollary 2.3 and Remark 2.4, our result actually allows “bi-level” corruptions/outliers that need to satisfy a much weaker requirement than this: the large-outliers have magnitude that is “large enough”, while the rest are such that the squared sum of their magnitudes is “small enough”. The threshold for both “large enough” and “small enough” decreases with each subspace update step.

II-C3 Order-wise improvement in allowed upper bound on maximum number of outliers per row

As pointed out in [6], solutions for standard RPCA (that only assume incoherence of left and right singular vectors of $\bm{L}$ and nothing else, i.e., no outlier support model) cannot tolerate666The reason is this: let $\rho_{\mathrm{row}}=\text{\small{max-outlier-frac-row}}$, one can construct a matrix $\bm{X}$ with $\rho_{\mathrm{row}}$ outliers in some rows that has rank equal to $1/\rho_{\mathrm{row}}$. A simple way to do this would be to let the support and nonzero entries of $\bm{X}$ be constant for $\rho_{\mathrm{row}}d$ columns before letting either of them change. Then the rank of $\bm{X}$ will be $d/(\rho_{\mathrm{row}}d)=1/\rho_{\mathrm{row}}$. If $1/\rho_{\mathrm{row}}<r_{\scriptscriptstyle{L}}=\operatorname{rank}(\bm{L})$, $\bm{X}$ will wrongly get classified as the low-rank component. This is why we need $\rho_{\mathrm{row}}=\text{\small{max-outlier-frac-row}}<1/r_{\scriptscriptstyle{L}}$. A similar argument can be used for max-outlier-frac-col. a bound on maximum outlier fractions in any row or any column that is larger than $1/r_{\scriptscriptstyle{L}}$. However observe that simple-ReProCS can tolerate $\text{\small{max-outlier-frac-row}}^{\alpha}\in O(1)$ (this assumes $f$ is a constant). This is a significant improvement over all existing RPCA results with important practical implications for video analytics. This is possible is because s-ReProCS uses extra assumptions, we explain their next.

II-C4 The need for extra assumptions

s-ReProCS recovers the sparse outliers $\bm{x}_{t}$ first and then the true data $\bm{\ell}_{t}$, and does this at each time $t$. Let $\bm{\Psi}:=\bm{I}-\hat{\bm{P}}_{(t-1)}\hat{\bm{P}}_{(t-1)}{}^{\prime}$. When recovering $\bm{x}_{t}$, it exploits two facts: (a) the subspace of $\bm{\ell}_{t}$, $\bm{P}_{(t)}$, satisfies the denseness/incoherence property, and (b) “good” knowledge of the subspace of $\bm{\ell}_{t}$ (either from initialization or from the previous subspace’s estimate and slow subspace change) is available. Using these two facts one can show that $\bm{\Psi}$ satisfies the RIP property, and that the “noise” seen by the compressive sensing step, $\bm{b}_{t}:=\bm{\Psi}(\bm{\ell}_{t}+\bm{v}_{t})$, is small. This, along with a guarantee for CS, helps ensure that the error in recovering $\bm{x}_{t}$ is upper bounded777Since the individual vector $\bm{b}_{t}$ does not have any structure that can be exploited, the error in recovering $\bm{x}_{t}$ cannot be made lower than this. However the $\bm{b}_{t}$’s arranged into a matrix do form a low-rank matrix whose approximate rank can be shown to be one (under our current subspace change model). If we try to exploit this structure we end up with the modified-PCP approach studied in earlier work [13]. This needs the uniform random support assumption [13]. by $C\|\bm{b}_{t}\|$. This, in turn, means that, to correctly recover the support of $\bm{x}_{t}$, the minimum large-outlier magnitude needs to be larger than $C\|\bm{b}_{t}\|$. This is where the $x_{\min}$ or $x_{\min,t}$ lower bound comes from888If there were a way to bound the element-wise error of the CS step (instead of the $l_{2}$ norm of the error), we could relax the $x_{\min}$ lower bound significantly. It is not clear if this is possible though..

Correct outlier support recovery is needed to ensure that the subspace estimate can be improved with each subspace update step. In particular, it helps ensure that the error vectors $\bm{e}_{t}:=\bm{x}_{t}-\bm{\hat{x}}_{t}$ in a given subspace update interval are mutually independent when conditioned on the $\bm{y}_{t}$’s from all past intervals. This fact also relies on the mutual independence assumption on the $\bm{a}_{t}$’s. Moreover, mutual independence, along with the element-wise boundedness and identical covariances assumption, on the $\bm{a}_{t}$’s helps ensure that we can use matrix Bernstein [31] and Vershynin’s sub-Gaussian result (bounds singular values of matrices with independent sub-Gaussian rows) [32] for obtaining the desired concentration bounds on the subspace recovery error in each step. As explained below, the above assumptions on $\bm{a}_{t}$ replace the right incoherence assumption. Finally, because ReProCS is updating the subspace using just the past $\alpha$ estimates of $\hat{\bm{\ell}}_{t}$’s, in order to show that each such step improves the estimate we need to bound $\text{\small{max-outlier-frac-row}}^{\alpha}$ (instead of just max-outlier-frac-row).

II-C5 Time and Memory Complexity

Observe that the s-ReProCS algorithm needs memory of order $n\alpha$ in online mode and $Kn\alpha$ in offline mode. Assuming $\alpha=\alpha_{*}$, even in offline mode, its memory complexity is near-optimal and equal to $O(nr\log n\log(1/\tilde{\varepsilon}))$. Also, observe that the time complexity of s-ReProCS is $O(ndr\log(1/\tilde{\varepsilon}))$. We explain this in Appendix D. These claims assume that $f$ is constant with $n,r$. If the dependence on $f$ is included both will be multiplied by $f^{2}$.

II-C6 Subspace and outlier assumptions’ tradeoff

When there are fewer outliers in the data or when outliers are easy to detect, one would expect to need weaker assumptions on the true data subspace or its rate of change. This is indeed true. For the original RPCA results, this is encoded in the condition $\max(\text{\small{max-outlier-frac-row}},\text{\small{max-outlier-frac-col}})\leq c/(\mu r_{\scriptscriptstyle{L}})$ where $\mu$ quantifies not-denseness of both left and right singular vectors. From Theorem 2.2, this is also how max-outlier-frac-col, $\mu$ (not-denseness of only left singular vectors) and $r_{\scriptscriptstyle{L}}$ are related for dynamic RPCA. On the other hand, for our result, $\text{\small{max-outlier-frac-row}}^{\alpha}$ and the lower bound on $x_{\min}$ govern the allowed rate of subspace change. The latter relation is easily evident from the bound on $\Delta$. If $x_{\min}$ is larger (outliers are large magnitude and hence easy to detect), a larger $\Delta$ can be tolerated. The relation of max-outlier-frac-row to rate of change is not evident from the way the guarantee is stated in Theorem 2.2. The reason is we have assumed $\text{\small{max-outlier-frac-row}}^{\alpha}\leq\rho_{\mathrm{row}}=0.01/f^{2}$ and used that to get a simple expression for $K$. If we did not do this, we would need $K$ to satisfy

[TABLE]

With this, $K$ needs to be $K=\lceil\frac{1}{-\log(c_{2}f\sqrt{\rho_{\mathrm{row}}})}\log(\frac{c_{1}\Delta}{0.8\tilde{\varepsilon}})\rceil$. Recall that we need $t_{j+1}-t_{j}\geq(K+3)\alpha$. Thus, a smaller $\rho_{\mathrm{row}}$ (smaller $\text{\small{max-outlier-frac-row}}^{\alpha}$) means one of two things: either a larger $\Delta$ (more change at each subspace change time) can be tolerated while keeping $K$, and hence the lower bound on the delay between change times, the same; or, for $\Delta$ fixed, a smaller lower bound is needed on the delay between change time. The above can be understood by carefully checking the proof999The multiplier 0.4 of $q_{\mathrm{rot}}$ in its first claim is obtained by setting $\rho_{\mathrm{row}}=0.01/f^{2}$. If we do not do this, 0.4 will get replaced by $c_{2}f\sqrt{\rho_{\mathrm{row}}}$. of Theorem 5.7.

II-C7 The need for detecting subspace change

As pointed out by an anonymous reviewer, it may not be clear to a reader why we need to explicitly detect subspace change (instead of just always doing subspace update at regularly spaced intervals). The change detection is needed for two key reasons. First, in the projection-SVD step for subspace update, we use $\hat{\bm{P}}_{j-1}$ as the best estimate of the previous subspace. We let $\hat{\bm{P}}_{j-1}$ be the final $\tilde{\varepsilon}$-accurate subspace estimate obtained after $K$ projection-SVD steps and then one subspace re-estimation step. To know when the $K$ updates are over, we need to know when the first update of the new subspace occurred, or in other words, we need an estimate of when the subspace change occurred. Second, in the current algorithm, because we detect change, we can choose to use the $\hat{\bm{\ell}}_{t}$’s from the next $\alpha$-frame interval, i.e. $[\hat{\bm{\ell}}_{{\hat{t}}_{j}+1},\hat{\bm{\ell}}_{{\hat{t}}_{j}+2},\dots,\hat{\bm{\ell}}_{{\hat{t}}_{j}+\alpha}]$, for the first subspace update. This ensures that the $\hat{\bm{\ell}}_{t}$’s from the interval that contains $t_{j}$ (some of the $\bm{\ell}_{t}$’s in this interval come from $\bm{P}_{j-1}$ while others come from $\bm{P}_{j}$) is never used further in any subspace update. The is essential because, if these are used, one will get an incorrect subspace estimate (something in between $\bm{P}_{j-1}$ and $\bm{P}_{j}$) and one whose error cannot easily be bounded. If subspace change is never detected, this cannot be ensured.

II-C8 Algorithm parameters

Observe from Theorem 2.2 that we need knowledge of only 4 model parameters - $r$, $\lambda^{+}$, $\lambda^{-}$ and $x_{\min}$ - to set our algorithm parameters. The initial dataset used for estimating $\hat{\bm{P}}_{0}$ (using PCP/AltProj) can also be used to get an accurate estimate of $r$, $\lambda^{-}$ and $\lambda^{+}$ using standard techniques (maximum likelihood applied to the AltProj estimate of $[\bm{\ell}_{1},\bm{\ell}_{2},\dots,\bm{\ell}_{t_{\mathrm{train}}}]$). Thus one really only needs to set $x_{\min}$. If continuity over time is assumed, a simple heuristic is to let it be time-varying and use $\min_{i\in\hat{\mathcal{T}}_{t-1}}|(\bm{\hat{x}}_{t-1})_{i}|$ as its estimate at time $t$. This approach in fact allows us to estimate $x_{\min,t}$ and thus allows for larger unstructured noise, $\bm{v}_{t}$, levels as allowed by Corollary 2.3.

The most interesting point for practice though is that of Remark 2.4. It indicates that when a subspace change is detected but not estimated, starting at the previous $2\alpha$ frames, one should use a larger value of the support estimation threshold $\omega_{supp}$. After each subspace update step, $\omega_{supp}$ should be decreased roughly exponentially.

II-C9 Dependence on $f$

Observe that $f$ appears in our guarantee in the bound on $\text{\small{max-outlier-frac-row}}^{\alpha}$ and in the expression for $\alpha$. The $\text{\small{max-outlier-frac-row}}^{\alpha}$ bound is stated that way only for simplicity. Actually, for all time instants except the $\alpha$-length period when the subspace re-estimation (for deletion) step is run, we only need $\text{\small{max-outlier-frac-row}}^{\alpha}\leq 0.01$. We need the tighter bound $\text{\small{max-outlier-frac-row}}^{\alpha}\leq\rho_{\mathrm{row}}=0.01/f^{2}$ only for the simple SVD based subspace re-estimation (deletion) step to work (i.e., only for $t\in[{\hat{t}}_{j}+K\alpha,{\hat{t}}_{j}+K\alpha+\alpha)$). Thus, if offline ReProCS were being used to solve the standard RPCA type problem (where $r_{\scriptscriptstyle{L}}$ is nicely bounded), one could choose to never run the subspace deletion step. This will mean that the resulting algorithm (s-ReProCS-no-delete) will need $\text{\small{max-outlier-frac-col}}<c/\mu r_{\scriptscriptstyle{L}}$, but then $\text{\small{max-outlier-frac-row}}<0.01$ will suffice (the bound would not depend on $1/f^{2}$). The $\alpha$ expression governs required delay between subspace change times, tracking delay, and time and memory complexity. If the deletion step is removed, the dependence of $\alpha$ on $f$ will not disappear, but will weaken (it will linearly depend on $f$ not on $f^{2}$).

We now try to relate $f$ to the condition number of $\bm{L}$. Observe that $f$ is the condition number of $\mathbb{E}[\bm{L}_{j}\bm{L}_{j}{}^{\prime}]$ for any $j$. The condition number of the entire matrix $\bm{L}$ can be much larger when slow subspace change holds ($\Delta$ is small). To see this, let $\kappa^{2}$ denote the condition number of $\mathbb{E}[\bm{L}\bm{L}^{\prime}]$, so that, whp, $\kappa$ is approximately the condition number of $\bm{L}$. It is not hard to to see that101010To understand this, suppose that there is only one subspace change and suppose that the intervals are equal, i.e., $d-t_{1}=t_{1}-t_{0}$. Then, $\mathbb{E}[\bm{L}\bm{L}^{\prime}]/d=\bm{P}_{0,{\mathrm{fix}}}\bm{\Lambda}_{\mathrm{fix}}\bm{P}_{0,{\mathrm{fix}}}{}^{\prime}+[\bm{P}_{0,{\mathrm{ch}}}\ \bm{P}_{1,\mathrm{new}}]\bm{B}[\bm{P}_{0,{\mathrm{ch}}}{}^{\prime}\ \bm{P}_{1,\mathrm{new}}{}^{\prime}]^{\prime}$ where \bm{B}=\lambda_{\mathrm{ch}}\left[\begin{array}[]{cc}(0.5+0.5\cos^{2}\theta_{1})&\ -0.5\sin\theta_{1}\cos\theta_{1}\\ -0.5\sin\theta_{1}\cos\theta_{1}&\ 0.5\sin^{2}\theta_{1}\end{array}\right]. The maximum eigenvalue of $\mathbb{E}[\bm{L}\bm{L}^{\prime}]/d$ is $\lambda^{+}$. Its minimum eigenvalue is the minimum eigenvalue of $\bm{B}$ which can be computed as $(1-\cos\theta_{1})\lambda_{\mathrm{ch}}$. In the worst case $\lambda_{\mathrm{ch}}=\lambda^{-}$. When the intervals are not equal, this gets replaced by $(1-\sqrt{1-2c\sin^{2}\theta_{1}})\lambda^{-}$ for a $c<1$. This is at most ${(1-\sqrt{1-2c\Delta^{2})}}$. , in the worst case (if $\lambda^{-}=\lambda_{\mathrm{ch}}$), $\kappa^{2}=\frac{f}{1-\sqrt{1-2c\Delta^{2}}}\approx C\frac{f}{\Delta^{2}}$ when $\Delta$ is small. Thus, if $\Delta$ is small, $\kappa\approx C\sqrt{f}/\Delta$ can be much larger than $f$. The guarantees of many of the RPCA solutions such as RPCA-GD [8] depend on $\kappa$.

II-C10 Relating our assumptions to right incoherence of $\bm{L}_{j}:=\bm{L}_{[t_{j},t_{j+1})}$ [4]

We repeat this discussion from [11]. From our assumptions, $\bm{L}_{j}=\bm{P}_{j}\bm{A}_{j}$ with $\bm{A}_{j}:=[\bm{a}_{t_{j}},\bm{a}_{t_{j}+1},\dots\bm{a}_{t_{j+1}-1}]$, the columns of $\bm{A}_{j}$ are zero mean, mutually independent, have identical covariance $\bm{\Lambda}$, $\bm{\Lambda}$ is diagonal, and are element-wise bounded as specified by Theorem 2.2. Let $d_{j}:=t_{j+1}-t_{j}$. Define a diagonal matrix $\Sigma$ with $(i,i)$-th entry $\sigma_{i}$ and with $\sigma_{i}^{2}:=\sum_{t}(a_{t})_{i}^{2}/d_{j}$. Define a $d_{j}\times r$ matrix $\tilde{\bm{V}}$ with the $t$-th entry of the $i$-th column being $(\tilde{\bm{v}}_{i})_{t}:=(\bm{a}_{t})_{i}/(\sigma_{i}\sqrt{d_{j}})$. Then, $\bm{L}_{j}=\bm{P}_{j}\Sigma\tilde{\bm{V}}^{\prime}$ and each column of $\tilde{\bm{V}}$ is unit 2-norm. Also, from the bounded-ness assumption, $(\tilde{\bm{v}}_{i})_{t}^{2}\leq\eta\frac{\lambda_{i}}{\sigma_{i}^{2}}\cdot\frac{1}{d_{j}}$ where $\eta$ is a numerical constant. Observe that $\bm{P}_{j}\Sigma\tilde{\bm{V}}^{\prime}$ is not exactly the SVD of $\bm{L}_{j}$ since the columns of $\tilde{\bm{V}}$ are not necessarily exactly mutually orthogonal. However, if $d_{j}$ is large enough, using the assumptions on $\bm{a}_{t}$, one can argue using any law of large numbers’ result (e.g., Hoeffding inequality), that (i) the columns of $\tilde{\bm{V}}$ are approximately mutually orthogonal whp, and (ii) $\sigma_{i}^{2}\geq 0.99\lambda_{i}$ whp. Thus, our assumptions imply that, whp, $\tilde{\bm{V}}$ is a basis matrix and $(\tilde{\bm{v}}_{i})_{t}^{2}\leq C/{d_{j}}$.

With the above, one can interpret $\tilde{\bm{V}}$ as an “approximation” to the right singular vectors of $\bm{L}_{j}$ and then the above bound on $(\tilde{\bm{v}}_{i})_{t}^{2}$ is the same as the right incoherence condition assumed by [4]. It is slightly stronger than what is assumed by [2, 6] and others (these do not require a bound on each entry but on each row, they require that the squared norm of each row of the matrix of right singular vectors be bounded by $Cr/d_{j}$).

Ideally we would like to work with the exact SVD of $\bm{L}_{j}$, however this is much harder to analyze using our statistical assumptions on the $\bm{a}_{t}$’s. To see this, suppose $\bm{A}_{j}\overset{\mathrm{SVD}}{=}\bm{U}\Sigma\bm{V}^{\prime}$, then $\bm{L}_{j}\overset{\mathrm{SVD}}{=}(\bm{P}_{j}\bm{U})\Sigma\bm{V}^{\prime}$ is the exact SVD of $\bm{L}_{j}$. Here $\bm{U}$ is an $r\times r$ orthonormal matrix. Now it is not clear how to relate the element-wise bounded-ness assumption on $\bm{a}_{t}$’s to an assumption on entries of $\bm{V}$, since now there is no easy expression for each entry of $\bm{V}$ or of the entries of $\Sigma$ in terms of $\bm{a}_{t}$ ($\bm{U}$ is an unknown matrix that can have all nonzero entries in general).

II-C11 Limitations of our guarantees

s-ReProCS needs a few extra assumptions beyond slow subspace change and what static RPCA solutions need: (i) instead of a bound on outlier fractions per row of the entire data matrix (which is what standard RPCA methods assume), it needs such a bound for every sub-matrix of $\alpha$ consecutive columns; (ii) it makes statistical assumptions on the principal subspace coefficients $\bm{a}_{t}$ (with mutual independence being the strongest requirement); (iii) it needs to lower bound $x_{\min}$; and (iv) it uses $\tilde{\varepsilon}$ to denote both the initial subspace error as well as the final recovery error achieved after a subspace update is complete. Here (i) is needed because ReProCS is an online algorithm that uses $\alpha$ frames at a time to update the subspace, and one needs to show that each update step provides an improved estimate compared to the previous one. However, since $\alpha$ is large enough, requiring a bound $\text{\small{max-outlier-frac-row}}^{\alpha}$ is not too much stronger than requiring the same bound on the outlier fractions per row of the entire $n\times d$ matrix $\bm{Y}$. In fact, if we compare the various RPCA solutions with storage complexity fixed at $O(n\alpha)=O(nr\log n)$, i.e., if we implement the various static RPCA solutions for every new batch of $\alpha$ frames of data, then, the static RPCA solutions will also need to bound $\text{\small{max-outlier-frac-row}}^{\alpha}$ defined in (8). As discussed earlier, these will require a much tighter bound of $c/r$ though. (iv) is assumed for simplicity. What we can actually prove is something slightly stronger: if the initial error is $\tilde{\varepsilon}$, and if $\epsilon=c\tilde{\varepsilon}$ for a constant $c$ which may be less than one, then, without any changes, we can guarantee the final subspace error to be below such an $\epsilon$. More generally, as long as the initial error $\tilde{\varepsilon}\leq\Delta$, it is possible to achieve final error $\epsilon$ for any $\epsilon>0$ if we assume that $t_{1}-t_{\mathrm{train}}>K\alpha$, assume a slightly larger lower bound on $x_{\min}$, and if we modify our initialization procedure to use the approach of follow-up work [11].

Limitations (ii) and (iii) are artifacts of our proof techniques. The mutual independence can be replaced by an autoregressive model on the $\bm{a}_{t}$’s by borrowing similar ideas from [7]. The mutual independence and zero mean assumption on the $\bm{a}_{t}$’s is valid for the video analytics’ application if we let $\bm{\ell}_{t}$ be the mean-subtracted background image at time $t$. Then, $\bm{\ell}_{t}$ models independent zero-mean background image variations about a fixed mean image, e.g., variations due to lighting variations or due to moving curtains; see Fig. 3. This type of mean subtraction (with an estimate of the mean background image computed from training data) is commonly done in practice in many practical applications where PCA is used; it is also done in our video experiments shown later. (iii) is needed because our proof first tries to show exact outlier support recovery by solving a CS problem to recover the outliers from the projected measurements, followed by thresholding. It should be possible to relax this by relaxing the exact support recovery requirement which, in turn, will require other significant changes. For example, it may be possible to do this if one is able to do a deterministic analysis. It may be possible to also completely eliminate it if we replace the CS step by thresholding with carefully decreasing thresholds in each iteration (borrow the idea of AltProj); however, we may then require the same tight bound on max-outlier-frac-row that AltProj needs. By borrowing the stagewise idea of AltProj, it may also be possible to remove all dependence on $f$.

III Discussion of Related Work

III-1 Limitations of earlier ReProCS-based guarantees [5, 12, 7]

In [5], we introduced the ReProCS idea and proved a partial guarantee for it. We call it a partial guarantee because it needed to assume something about the intermediate subspace estimates returned by the algorithm. However, this work is important because it developed a nice framework for proving guarantees for dynamic RPCA solutions. Both our later complete guarantees [12, 7] as well as the current result build on this framework.

The current work is a significant improvement over the complete guarantees obtained in [12, 7] for two other ReProCS-based algorithms for three reasons. (i) The earlier works needed very specific assumptions on how the outlier support could change (needed an outlier support model inspired by video moving objects). Our result removes such a requirement and instead only needs a bound on the fraction of outliers per column of the data matrix and on the fraction per row of an $\alpha$-consecutive-column sub-matrix of the data matrix (for $\alpha$ large enough). (ii) The subspace change model assumed in these earlier papers can be interpreted as the current model (given in Sec. II) with $\theta_{j}=90^{\circ}$ or equivalently with $\Delta=1$. This is an unrealistic model for slow subspace change, e.g., in 3D, it implies that the subspace changes from the x-y plane to the y-z plane. Instead, our current model allows changes from x-y plane to a slightly tilted x-y plane as shown in Fig. 1. This modification is more realistic and it allows us to replace the upper bound on $\lambda_{\mathrm{ch}}$ required by the earlier results by a similar bound on $\lambda_{\mathrm{ch}}\Delta^{2}$ (see assumption 1b of Theorem 2.2). Since $\Delta$ quantifies rate of subspace change, this new requirement is much weaker. It can be satisfied by assuming that $\Delta$ is small, without making any assumption on $\lambda_{\mathrm{ch}}$. (iii) The required minimum delay between subspace change times in the earlier results depended on $1/\epsilon^{2}$ where $\epsilon$ is the desired final subspace error after a subspace update is complete. This is a strong requirement. Our current result removes this unnecessarily strong dependence. The delay now only depends on $(-\log\epsilon)$ which makes it much smaller. It also implies that the memory complexity of simple-ReProCS is near-optimal. (iv) Unlike [12, 7], we analyze a simple ReProCS-based algorithm that ensures that the estimated subspace dimension is bounded by $(r+1)$, without needing the complicated cluster-SVD algorithm. This is why our guarantee allows outlier fractions per column to be below $c/r$. The work of [12] needed this to be below $c/r_{\scriptscriptstyle{L}}$ while [7] needed an extra assumption (clustered eigenvalues). For long data sequences, $c/r$ can be much larger than $c/r_{\scriptscriptstyle{L}}$. We provide a detailed comparison of these assumptions in Table II.

III-2 Complete guarantees for other dynamic RPCA or RPCA solutions

Another approach that solves dynamic RPCA, but in a piecewise batch fashion, is modified-PCP (mod-PCP) [13]. The guarantee for mod-PCP was proved using ideas borrowed from [2] for PCP. Thus, like [2], it also needs uniformly randomly generated support sets which is an unrealistic requirement. For the video application, this requires that foreground objects are single pixel wide and move around the entire image completely randomly over time. This is highly impractical. In Table I, we provide a comparison of our current guarantees for simple-ReProCS and its offline version with those for original-ReProCS [5, 12, 7], modified-PCP [13], as well as with those for solutions for standard RPCA - [2] (referred to as PCP(C)), [4] (PCP(H), this strictly improves upon [3]), AltProj [6], RPCA via gradient descent (GD) [8] and nearly-optimal robust matrix completion (NO-RMC) [18]. The table also contains a speed and memory complexity comparison. Offline s-ReProCS can be interpreted as a solution for standard RPCA. From Table I, it is clear that for data that satisfies slow subspace change and the assumption that outlier magnitudes are either large or very small, and that is such that its first $t_{\mathrm{train}}$ frames, $\bm{Y}_{[1,t_{\mathrm{train}}]}$, satisfy AltProj (or PCP) assumptions, s-ReProCS and offline s-ReProCS have the following advantages over other methods.

For the data matrix after $t_{\mathrm{train}}$, i.e., for $\bm{Y}_{[t_{\mathrm{train}}+1,d]}$, ReProCS needs the weakest bound on $\text{\small{max-outlier-frac-row}}^{\alpha}$ without requiring uniformly randomly generated outlier support sets. This is comparable to the bound needed by PCP(C) or mod-PCP but both assume uniform random outlier supports which is a very strong requirement. 2. 2.

The memory complexity of s-ReProCS is significantly better than that of all other published methods for RPCA that provably work, and is nearly optimal. 3. 3.

Both in terms of time complexity order (Table I) and experimentally (see Sec. VI), s-ReProCS and its offline counterpart are among the fastest, while having the best, or nearly the best, performance experimentally as well. Order-wise, only NO-RMC [18] is faster than s-ReProCS. However, NO-RMC needs the data matrix to be nearly square, i.e., it needs $c_{1}n\geq d\geq c_{2}n$. This is a very strong requirement that often does not hold: for the video application it requires that the number of video frames $d$ be roughly as large as $n$ (number of pixels in one image frame). The reason is that NO-RMC deliberately under-samples the data matrix $\bm{Y}$ by randomly throwing away some of its entries and using only the rest even when all are available. In other words, it always solves the robust matrix completion (RPCA with missing entries) problem and this is what results in a significant speed-up, but this is also why it needs $d\approx n$. 4. 4.

s-ReProCS can automatically detect subspace change and then also track it with a short delay, while the other approaches (except original-ReProCS) cannot. Notice that s-ReProCS also needs a weaker upper bound of $c/r$ on max-outlier-frac-col while the batch techniques (PCP, AltProj, GD, NO-RMC) applied to the entire matrix $\bm{L}$ need this to be below $c/r_{\scriptscriptstyle{L}}$. Of course, if the batch techniques are applied on pieces of data $\bm{Y}_{j}:=\bm{Y}_{[t_{j},t_{j+1})}$, they also need the same looser bound of $c/r$ on max-outlier-frac-col and their memory complexity improves too. However, the batch methods do not have a way to estimate the change times $t_{j}$, while s-ReProCS does. Moreover, since the other methods (except modified-PCP) do not have a way to use the previous subspace information, if the pieces chosen are are too small, e.g., if the methods are applied on $\alpha$-frames at a time, their performance is much worse then when the entire dataset is used jointly.

III-3 Comparison with follow-up work on ReProCS-NORST [11]

As compared to ReProCS-NORST, which is the algorithm studied in our follow-up work [11] (which allows all $r$ directions of the subspace to change at each $t_{j}$), simple-ReProCS has three advantages: (i) it is faster, (ii) it needs a weaker lower bound on $x_{\min}$ (its required lower bound essentially does not depend on $r$ if $\tilde{\varepsilon}$ is very small), and (iii) if it is used to solve the standard RPCA problem (estimate span of columns of the entire matrix $\bm{L}$), we can eliminate the $r$-SVD based subspace re-estimation (deletion) step. With this change, (a) the required upper bound on $\text{\small{max-outlier-frac-row}}^{\alpha}$ for s-ReProCS does not depend on the condition number $f$ (just $\text{\small{max-outlier-frac-row}}^{\alpha}\leq 0.01$ suffices), and (b) its time complexity improves by a factor of $r$ (if the initialization step is ignored) compared to ReProCS-NORST. Of course it will mean s-ReProCS-no-delete will need $\text{\small{max-outlier-frac-col}}<c/(r+J)$ which is slightly stronger.

Since ReProCS-NORST allows all $r$ directions of the subspace to change, it also has many advantages over s-ReProCS: its subspace tracking delay is near-optimal, and it allows for a weaker initialization assumption.

IV Why s-ReProCS works: main ideas of our proof

In this section we explain the main ideas of our proof, first for the $t_{j}$ known case, and then explain why the subspace change detection step works.

IV-A Why s-ReProCS with $t_{j}$ known works

To understand things simply, first assume that ${\hat{t}}_{j}=t_{j}$, i.e., the subspace change times are known. Consider Algorithm 1. At each time $t$ this consists of three steps - projected Compressive Sensing (CS) to estimate $\bm{x}_{t}$, estimating $\bm{\ell}_{t}$ by subtraction, and subspace update. Consider projected CS. This is analyzed in Lemma 5.15. At time $t$, suppose that we have access to $\hat{\bm{P}}_{(t-1)}$ which is a good estimate of the previous subspace, $\operatorname{span}(\bm{P}_{(t-1)})$. Because of slow subspace change, this is also a good estimate of $\operatorname{span}(\bm{P}_{(t)})$. Its first step projects $\bm{y}_{t}$ orthogonal to $\hat{\bm{P}}_{(t-1)}$ to get $\tilde{\bm{y}}_{t}$. Recall that $\tilde{\bm{y}}_{t}=\bm{\Psi}\bm{x}_{t}+\bm{b}_{t}$ where $\bm{b}_{t}:=\bm{\Psi}(\bm{\ell}_{t}+\bm{v}_{t})$ is small. Using the incoherence (denseness) assumption and $\operatorname{span}(\bm{P}_{(t-1)})$ being a good estimate of $\operatorname{span}(\bm{P}_{(t)})$, it can be argued that the restricted isometry constant (RIC) [10] of $\bm{\Psi}:=\bm{I}-\hat{\bm{P}}_{(t-1)}\hat{\bm{P}}_{(t-1)}{}^{\prime}$ will be small. Using [10, Theorem 1.2], this, along with $\|\bm{b}_{t}\|$ being small, ensures that $l_{1}$ minimization will produce an accurate estimate, $\bm{\hat{x}}_{t,cs}$, of $\bm{x}_{t}$. The support estimation step with a carefully chosen threshold, $\omega_{supp}=x_{\min}/2$, and a lower bound on $x_{\min}$ then ensures exact support recovery, i.e., $\hat{\mathcal{T}}_{t}=\mathcal{T}_{t}$. With this, the LS step output, $\bm{\hat{x}}_{t}$, satisfies $\bm{\hat{x}}_{t}=\bm{x}_{t}+\bm{e}_{t}$ with

[TABLE]

and with $\|\bm{e}_{t}\|$ being small. Computing $\hat{\bm{\ell}}_{t}:=\bm{y}_{t}-\bm{\hat{x}}_{t}$, then also gives a good estimate of $\bm{\ell}_{t}$ that satisfies $\hat{\bm{\ell}}_{t}=\bm{\ell}_{t}+\bm{v}_{t}-\bm{e}_{t}$ with $\bm{e}_{t}$ as above.

The subspace update step uses $\hat{\bm{\ell}}_{t}$’s to update the subspace. Since $\bm{e}_{t}$ satisfies (10), $\bm{e}_{t}$ depends on $\bm{\ell}_{t}$; thus the error/noise, $\bm{v}_{t}-\bm{e}_{t}$, in the “observed data” $\hat{\bm{\ell}}_{t}$ used for the subspace update step depends on the true data $\bm{\ell}_{t}$. Because of this, the subspace update does not involve a PCA or an incremental PCA problem in the traditionally studied setting (data and corrupting noise/error being independent or uncorrelated). It is, in fact, an instance of PCA when the noise/error, $\bm{v}_{t}-\bm{e}_{t}$, in the observed data $\hat{\bm{\ell}}_{t}$ depends on the true data $\bm{\ell}_{t}$. This problem was studied in [22, 23] where it was referred to as “correlated-PCA” or “PCA in data-dependent noise”. Using this terminology, our subspace update problem (estimating $\bm{P}_{j}$ using $\hat{\bm{P}}_{j-1}$) is a problem of PCA in data-dependent noise with partial subspace knowledge. To simplify our analysis, we first study this more general problem and obtain a guarantee for it in Theorem 5.7 in Sec. V-A. This theorem along with Lemma 5.15 (that analyzes the projected-CS step discussed above) help obtain a guarantee for the $k$-th projection-SVD step in Lemma 5.16. The $k=1$ and $k>1$ cases are handled separately. The main assumption required for applying Theorem 5.7 holds because $\bm{e}_{t}$ is sparse with support $\mathcal{T}_{t}$ that changes enough ($\text{\small{max-outlier-frac-row}}^{\alpha}$ bound of Theorem 2.2 holds). The subspace deletion via simple SVD step of subspace update is studied in Lemma 5.17. This step solves a problem of PCA in data-dependent noise and so it directly uses the results from [23].

To understand the flow of the proof, consider the interval $[t_{j},t_{j+1})$. Assume that, before $t_{j}$, the previous subspace has been estimated with error $\tilde{\varepsilon}$, i.e., we have $\hat{\bm{P}}_{j-1}$ with $\mathrm{SE}(\hat{\bm{P}}_{j-1},\bm{P}_{j-1})\leq\tilde{\varepsilon}$. We explain below that this implies that, under the theorem’s assumptions, we get $\mathrm{SE}(\hat{\bm{P}}_{j},\bm{P}_{j})\leq\tilde{\varepsilon}$ before $t_{j+1}$. We remove the subscripts $j$ in some of this discussion. Define the interval $\mathcal{J}_{k}:=[t_{j}+(k-1)\alpha,t_{j}+k\alpha)$. Suppose also that $\bm{v}_{t}=0$.

Before the first projection-SVD step (which is done at $t=t_{j}+\alpha$), i.e., for $t\in\mathcal{J}_{1}$, we have no estimate of $\bm{P}_{\mathrm{new}}$, and hence only a crude estimate of $\bm{P}_{{\mathrm{rot}}}$. In particular, we can only get the bound $\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{\mathrm{rot}})=\mathrm{SE}(\hat{\bm{P}}_{j-1},\bm{P}_{\mathrm{rot}})\leq\tilde{\varepsilon}+|\sin\theta|$ for this interval.

• •

As a result, the bound on the “noise”, $\bm{b}_{t}$, seen by the projected-CS step is also the largest for this interval, we have $\|\bm{b}_{t}\|\leq C(\tilde{\varepsilon}\sqrt{r\lambda^{+}}+|\sin\theta|\sqrt{\lambda_{\mathrm{ch}}})$. Using the CS guarantee, followed by ensuring exact support recovery (as explained above), this implies that $\bm{e}_{t}$ satisfies (10) and that we get a similar bound on the final CS step error: $\|\bm{e}_{t}\|\leq C(\tilde{\varepsilon}\sqrt{r\lambda^{+}}+0.11|\sin\theta|\sqrt{\lambda_{\mathrm{ch}}})$. The factor of $0.11$ in the second term of this bound is obtained because, for this interval, $\bm{\Psi}=\bm{I}-\hat{\bm{P}}_{j-1}\hat{\bm{P}}_{j-1}{}^{\prime}$ and so $\bm{\Psi}\bm{P}_{\mathrm{new}}\approx\bm{P}_{\mathrm{new}}$ and $\bm{P}_{\mathrm{new}}$ is dense, see (13). Thus one can show that $\|\bm{I}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}\bm{P}_{\mathrm{new}}\|_{2}\leq 0.11$.

• •

This bound on $\bm{e}_{t}$, along with using the critical fact that $\bm{e}_{t}$ satisfies (10) (is sparse) and its support $\mathcal{T}_{t}$ changes enough (the $\text{\small{max-outlier-frac-row}}^{\alpha}$ bound of Theorem 2.2 holds), ensures that we get a better estimate of $\bm{P}_{\mathrm{rot}}$ after the first projection-SVD step. This is what allows us to apply Theorem 5.7. Using it we can show that $\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{\mathrm{rot}})=\mathrm{SE}([\hat{\bm{P}}_{j-1},\hat{\bm{P}}_{j,{\mathrm{rot}},1}],\bm{P}_{\mathrm{rot}})\leq 0.1\tilde{\varepsilon}+0.06|\sin\theta|$ for $t\in\mathcal{J}_{2}$. See proof of $k=1$ case of Lemma 5.16 and Fact 5.14. 2. 2.

Thus we have a much better estimate of $\bm{P}_{\mathrm{rot}}$ for $t\in\mathcal{J}_{2}$ than for $\mathcal{J}_{1}$. Because of this, $\|\bm{b}_{t}\|$ is smaller, and hence $\|\bm{e}_{t}\|$ is smaller for $t\in\mathcal{J}_{2}$. This, along with the sparsity and changing support, $\mathcal{T}_{t}$, of $\bm{e}_{t}$, ensures an even better estimate at the second projection-SVD step. We can show that $\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{\mathrm{rot}})=\mathrm{SE}([\hat{\bm{P}}_{j-1},\hat{\bm{P}}_{j,{\mathrm{rot}},2}],\bm{P}_{\mathrm{rot}})\leq 0.1\tilde{\varepsilon}+0.5\cdot 0.06|\sin\theta|$ for $t\in\mathcal{J}_{3}$. See proof of $k>1$ case of Lemma 5.16 and Fact 5.14. 3. 3.

Proceeding this way, we show that $\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{\mathrm{rot}})=\mathrm{SE}([\hat{\bm{P}}_{j-1},\hat{\bm{P}}_{j,{\mathrm{rot}},k}],\bm{P}_{\mathrm{rot}})\leq 0.1\tilde{\varepsilon}+0.5^{k-2}(0.06|\sin\theta|)$ after the $k$-th projection-SVD step. Picking $K$ appropriately, gives $\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{\mathrm{rot}})\leq\tilde{\varepsilon}$ after $K$ steps, i.e., at $t=t_{j}+K\alpha$. In all the above intervals, $\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{(t)})\leq\tilde{\varepsilon}+\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{\mathrm{rot}})$. Thus at $t=t_{j}+K\alpha$, $\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{(t)})\leq 2\tilde{\varepsilon}$. 4. 4.

At $t=t_{j}+K\alpha$, $\hat{\bm{P}}_{(t)}$ contains $(r+1)$ columns. The subspace re-estimation via simple SVD step re-estimates $\bm{P}_{j}$ in order to delete the deleted direction, $\bm{P}_{\mathrm{del}}$, from $\hat{\bm{P}}_{(t)}$. The output of this step is $\hat{\bm{P}}_{j}$ (the final estimate of $\operatorname{span}(\bm{P}_{j})$). Thus, at $t=t_{j}+K\alpha+\alpha$, $\hat{\bm{P}}_{(t)}=\hat{\bm{P}}_{j}$ and we can show that it satisfies $\mathrm{SE}(\hat{\bm{P}}_{(t)},\bm{P}_{(t)})=\mathrm{SE}(\hat{\bm{P}}_{j},\bm{P}_{j})\leq\tilde{\varepsilon}$. See Lemma 5.17. The re-estimation is done at this point because, for times $t$ in this interval, $\|\hat{\bm{\ell}}_{t}-\bm{\ell}_{t}\|=\|\bm{e}_{t}\|\leq 2.4\tilde{\varepsilon}\|\bm{\ell}_{t}\|$. For PCA in data-dependent noise, simple SVD needs $\alpha\geq(q/\epsilon)^{2}f^{2}(r\log n)$ where $q$ is the error/noise to signal ratio and $\epsilon$ is the final desired error level. For our problem, the “noise” is $\bm{e}_{t}$ and thus $q=2.4\tilde{\varepsilon}$ and $\epsilon=\tilde{\varepsilon}$. Since $q/\epsilon$ is a constant, $\alpha\geq\alpha_{*}=Cf^{2}r\log n$ suffices when simple SVD is applied at this time.

When $\bm{v}_{t}\neq 0$, almost all of the above discussion remains the same. The reason is this: in the main theorem, we assume $\|\bm{v}_{t}\|^{2}\leq 0.1r\tilde{\varepsilon}^{2}\lambda^{+}$ with $c<1$ and so even though we have to deal with $\bm{v}_{t}$ in $\bm{b}_{t}$ and $\bm{e}_{t}$ expressions, and in the $\alpha$ expression, the changes required are only to the numerical constants. In Corollary 2.3, we have carefully chosen the bound on $\|\bm{v}_{t}\|$ to equal the bound on $\|\bm{e}_{t}\|$ modulo constants. Thus, once again, only numerical constants change, everything else remains the same.

IV-B Why automatic subspace change detection and Automatic Simple-ReProCS works

The subspace change detection approach is summarized in Algorithm 4. This idea is motivated by a similar idea first used in our earlier works [12, 7]. The algorithm toggles between the “detect” phase and the “update” phase. It starts in the “detect” phase. If the $j$-th subspace change is detected at time $t$, we set ${\hat{t}}_{j}=t$. At this time, the algorithm enters the “update” (subspace update) phase. We then repeat the $K$ projection-SVD steps and the one subspace re-estimation via simple SVD step from Algorithm 1 with the following change: the $k$-th projection-SVD step is now done at $t={\hat{t}}_{j}+k\alpha-1$ (instead of at $t=t_{j}+k\alpha-1$) and the subspace re-estimation is done at $t={\hat{t}}_{j}+K\alpha+\alpha-1:={\hat{t}}_{j,fin}$. Thus, at $t={\hat{t}}_{j,fin}$, the subspace update is complete. At this time, the algorithm enters the “detect” phase again. To understand the change detection strategy, consider the $j$-th subspace change. Assume that the previous subspace $\bm{P}_{j-1}$ has been accurately estimated by $t={\hat{t}}_{j-1,fin}$ and that ${\hat{t}}_{j-1,fin}<t_{j}$. Let $\hat{\bm{P}}_{*}:=\hat{\bm{P}}_{j-1}$ denote this estimate. At this time, the algorithm enters the “detect” phase in order to detect the next ($j$-th) change. Let $\bm{B}_{t}:=(\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime})[\hat{\bm{\ell}}_{t-\alpha+1},\dots,\hat{\bm{\ell}}_{t}]$. For every $t={\hat{t}}_{j-1,fin}+u\alpha$, $u=1,2,\dots$, we detect change by checking if the maximum singular value of $\bm{B}_{t}$ is above a pre-set threshold, $\sqrt{\omega_{evals}\alpha}$, or not.

We claim that, whp, under Theorem 2.2 assumptions, this strategy has no false detects and correctly detects change within a delay of at most $2\alpha$ frames. The former is true because, for any $t$ for which $[t-\alpha+1,t]\subseteq[{\hat{t}}_{j-1,fin},t_{j})$, all singular values of the matrix $\bm{B}_{t}$ will be close to zero (will be of order $\sqrt{\tilde{\varepsilon}}$) and hence its maximum singular value will be below $\sqrt{\omega_{evals}\alpha}$. Thus, whp, ${\hat{t}}_{j}\geq t_{j}$. To understand why the change is correctly detected within $2\alpha$ frames, first consider $t={\hat{t}}_{j-1,fin}+\left\lceil\frac{t_{j}-{\hat{t}}_{j-1,fin}}{\alpha}\right\rceil\alpha:=t_{j,*}$. Since we assumed that ${\hat{t}}_{j-1,fin}<t_{j}$ (the previous subspace update is complete before the next change), $t_{j}$ lie in the interval $[t_{j,*}-\alpha+1,t_{j,*}]$. Thus, not all of the $\bm{\ell}_{t}$’s in this interval will be generated from $\operatorname{span}(\bm{P}_{j})$. Thus, depending on where in the interval $t_{j}$ lies, the algorithm may or may not detect the change at this time. However, in the next interval, i.e., for $t\in[t_{j,*}+1,t_{j,*}+\alpha]$, all of the $\bm{\ell}_{t}$’s will be generated from $\operatorname{span}(\bm{P}_{j})$. We can prove that, whp, $\bm{B}_{t}$ for this time $t$ will have maximum singular value that is above the threshold. Thus, if the change is not detected at $t_{j,*}$, whp, it will get detected at $t_{j,*}+\alpha$. Hence one can show that, whp, either ${\hat{t}}_{j}=t_{j,*}$, or ${\hat{t}}_{j}=t_{j,*}+\alpha$, i.e., $t_{j}\leq{\hat{t}}_{j}\leq t_{j}+2\alpha$. To see the actual proof of these claims, please refer to Appendix A where we prove our main result without assuming $t_{j}$ known.

V Proving Theorem 2.2 with assuming ${\hat{t}}_{j}=t_{j}$

In Table IV, we summarize all the new symbols and terms used in our proof. This, along with Table III given earlier, should help follow the proof details without having to refer back to earlier sections. To give a simpler proof first, we prove Theorem 2.2 under the assumption that ${\hat{t}}_{j}=t_{j}$ below. The proof without this assumption is given in Appendix A. With assuming ${\hat{t}}_{j}=t_{j}$, we are studying Algorithm 1. Recall from Sec. IV that the subspace update step involves solving a problem of PCA in data-dependent noise when partial subspace knowledge is available. We provide a guarantee for this problem in Sec. V-A and use it in our proof in Sec. V-D.

For the entire proof we will use the equivalent subspace change model described in (5). Clearly $|\sin\theta_{j}|\leq\Delta$ by our assumption.

V-A PCA in data-dependent noise with partial subspace knowledge

V-A1 The Problem

We are given a set of $\alpha$ frames of observed data $\bm{y}_{t}:=\bm{\ell}_{t}+\bm{w}_{t}+\bm{z}_{t}$, with $\|\bm{z}_{t}\|^{2}\leq b_{z}^{2}$ and $\|\mathbb{E}[\bm{z}_{t}\bm{z}_{t}{}^{\prime}]\|\leq\lambda_{z}^{+}:=c_{1}b_{z}^{2}/r$ ; $\bm{w}_{t}=\bm{M}_{t}\bm{\ell}_{t}$, $\bm{\ell}_{t}=\bm{P}\bm{a}_{t}$ and with $\bm{P}$ satisfying

[TABLE]

$\bm{P}_{{\mathrm{fix}}}=(\bm{P}_{*}\bm{U}_{0})\bm{I}_{[1,r-1]}$, $\bm{P}_{{\mathrm{ch}}}=(\bm{P}_{*}\bm{U}_{0})\bm{I}_{r}$, and $\bm{U}_{0}$ is an $r\times r$ rotation matrix. Also, the $\bm{a}_{t}$’s are zero mean, mutually independent, element-wise bounded r.v.’s with identical and diagonal covariance $\bm{\Lambda}$, and independent of the matrices $\bm{M}_{t}$. The matrices $\bm{M}_{t}$ are unknown.

Let $\mathcal{J}^{\alpha}$ denote the $\alpha$-frame time interval for which the $\bm{y}_{t}$’s are available. We also have access to a partial subspace estimate $\hat{\bm{P}}_{*}$ that satisfies $\mathrm{SE}(\hat{\bm{P}}_{*},\bm{P}_{*})\leq\tilde{\varepsilon}$, and that is computed using data that is independent of the $\bm{\ell}_{t}$’s (and hence of the $\bm{y}_{t}$’s) for $t\in\mathcal{J}^{\alpha}$. The goal is to estimate $\operatorname{span}(\bm{P})$ using $\hat{\bm{P}}_{*}$ and the $\bm{y}_{t}$’s for $t\in\mathcal{J}^{\alpha}$.

V-A2 Projection-SVD / Projection-EVD

Let $\bm{\Phi}:=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}$. A natural way to estimate $\bm{P}$ is to first compute $\hat{\bm{P}}_{\mathrm{rot}}$ as the top eigenvector of

[TABLE]

and set $\hat{\bm{P}}=[\hat{\bm{P}}_{*},\hat{\bm{P}}_{\mathrm{rot}}]$. We refer to this strategy as “projection-EVD” or “projection-SVD”. In this paper, we are restricting ourselves to only one changed direction and hence we compute only the top eigenvector (or left singular vector) of $\bm{D}_{obs}$. In general if there were $r_{\mathrm{ch}}>1$ directions, we would compute all eigenvectors with eigenvalues above a threshold, see [5, 7].

V-A3 The Guarantee

We can prove the following about projection-SVD.

Theorem 5.7**.**

Consider the above setting for an $\alpha\geq\alpha_{0}$ where

[TABLE]

Assume that the $\bm{M}_{t}$’s can be decomposed as $\bm{M}_{t}=\bm{M}_{2,t}\bm{M}_{1,t}$ where $\bm{M}_{2,t}$ is such that $\left\|\bm{M}_{2,t}\right\|\leq 1$ but

[TABLE]

Let $q_{0}$ denote a bound on $\max_{t}\|\bm{M}_{1,t}\bm{P}_{*}\|$ and let $q_{\mathrm{rot}}$ denote a bound on $\max_{t}\|\bm{M}_{1,t}\bm{P}_{\mathrm{rot}}\|$, i.e., we have $\|\bm{M}_{1,t}\bm{P}_{*}\|\leq q_{0}\ \text{ and }\ \|\bm{M}_{1,t}\bm{P}_{\mathrm{rot}}\|\leq q_{\mathrm{rot}}$ for all $t\in\mathcal{J}^{\alpha}$. Assume that

[TABLE]

Define the event $\mathcal{E}_{*}:=\{\mathrm{SE}(\hat{\bm{P}}_{*},\bm{P}_{*})\leq\tilde{\varepsilon}\}$. The following hold.

*Conditioned on $\mathcal{E}_{*}$, w.p. at least $1-12n^{-12}$, *

$\mathrm{SE}(\hat{\bm{P}},\bm{P})\leq\tilde{\varepsilon}+\mathrm{SE}(\hat{\bm{P}},\bm{P}_{\mathrm{rot}})$* and*

[TABLE] 2. 2.

Conditioned on $\mathcal{E}_{*}$, w.p. at least $1-12n^{-12}$,

[TABLE]

For large $n,r$, $r\log n>r+\log n$. Thus the following simpler expression for $\alpha_{0}$ suffices: $\alpha\geq\alpha_{0}=C\eta^{2}f^{2}(r\log n)$.

Remark 5.8**.**

Theorem 5.7 holds even when $\mathcal{E}_{*}$ is replaced by $\mathcal{E}_{0}:=\mathcal{E}_{*}\cap\tilde{\mathcal{E}}(Z)$ where $\tilde{\mathcal{E}}(Z)$ is an event that depends on a r.v. $Z$ that is such that the pair $\{\hat{\bm{P}}_{*},Z\}$ is still independent of the $\bm{\ell}_{t}$’s (and hence of the $\bm{y}_{t}$’s) for $t\in\mathcal{J}^{\alpha}$.

Proof.

The proof follows using a careful application of the Davis-Kahan $\sin\theta$ theorem [33] followed by using matrix Bernstein [31] to bound the numerator terms in the $\sin\theta$ theorem bound and Vershynin’s sub-Gaussian result [32] to bound the extra terms in its denominator. While the overall approach is similar to that used by [23] for the basic correlated-PCA problem, this proof requires significantly more work. We give the proof in Appendix B. The most important idea in the proof is the use of Cauchy-Schwarz to show that the time-averaged projected-data - noise correlation and time-averaged noise power are both $\sqrt{b_{0}}$ times their instantaneous values. We explain this next. ∎

In the result above, the bounds assumed in (12) are not critical. They only help to get a simple expression for the subspace error bound. As will be evident from the proof, we can also get a (more complicated) guarantee without assuming (12), and with any value of $\rho_{\mathrm{row}}$. The main assumption needed by Theorem 5.7 is (11) on the data-dependency matrices $\bm{M}_{t}$. This is required because the noise $\bm{w}_{t}$ depends on the true data $\bm{\ell}_{t}$ and hence the instantaneous values of both the noise power and of the signal-noise correlation (even after being projected orthogonal to $\hat{\bm{P}}_{*}$) can be large if $\lambda_{\mathrm{ch}}$ is large. However, (11) helps ensure that the time-averaged noise power and the time-averaged projected-signal-noise correlation are much smaller. Using the definitions of $q_{0}$ and $q_{\mathrm{rot}}$, $\|\mathbb{E}[\bm{w}_{t}\bm{w}_{t}{}^{\prime}]\|\leq q_{0}^{2}\lambda^{+}+q_{\mathrm{rot}}^{2}\lambda_{\mathrm{ch}}:=c_{w}$ and $\|\mathbb{E}[\bm{\Phi}\bm{\ell}_{t}\bm{w}_{t}{}^{\prime}]\|\leq\tilde{\varepsilon}q_{0}\lambda^{+}+(\tilde{\varepsilon}+|\sin\theta|)q_{\mathrm{rot}}\lambda_{\mathrm{ch}}:=c_{wl}$. By appropriately applying Cauchy-Schwarz (Theorem E.26),

[TABLE]

[TABLE]

Since $b_{0}=0.01$, both bounds are $0.1$ times their instantaneous values. See proof of Lemma B.24 for their derivation.

For our problem, (11) holds because we can let $\bm{M}_{2,t}=\bm{I}_{\mathcal{T}_{t}}$ and $\bm{M}_{1,t}$ as the rest of the matrix multiplying $\bm{\ell}_{t}$ in (10). Then, using the bound on outlier fractions per row from Theorem 2.2, it is not hard to see that $\left\|\frac{1}{\alpha}\sum_{t}\bm{M}_{2,t}\bm{M}_{2,t}{}^{\prime}\right\|\leq\rho_{\mathrm{row}}$. We state this formally next in Lemma 5.9.

Lemma 5.9**.**

Assume that the $\text{\small{max-outlier-frac-row}}^{\alpha}$ bound of Theorem 2.2 holds. Then, for any $\alpha$-length interval $\mathcal{J}^{\alpha}\subseteq[t_{1},d]$,

[TABLE]

Proof of Lemma 5.9.

The proof is straightforward. Let $\bm{C}_{t}:=\bm{I}_{\mathcal{T}_{t}}\bm{I}_{\mathcal{T}_{t}}{}^{\prime}$. Then,

[TABLE]

where

[TABLE]

Since each $\bm{C}_{t}$ is diagonal, so is $\frac{1}{\alpha}\sum_{t\in\mathcal{J}^{\alpha}}\bm{C}_{t}$. The latter has diagonal entries given by

[TABLE]

Thus,

[TABLE]

where $\gamma(\mathcal{J}^{\alpha})$ defined in (7) is the outlier fraction per row of $\bm{X}_{\mathcal{J}^{\alpha}}$. From Theorem 2.2, this is bounded by $\rho_{\mathrm{row}}$. ∎

V-B Two simple lemmas from [5]

The following two lemmas taken from [5] will be used in proving the main lemmas that together imply Theorem 2.2 with ${\hat{t}}_{j}=t_{j}$.

Lemma 5.10**.**

[5, Lemma 2.10]** Suppose that $\bm{P}$, $\hat{\bm{P}}$ and $\bm{Q}$ are three basis matrices. Also, $\bm{P}$ and $\hat{\bm{P}}$ are of the same size, $\bm{Q}^{\prime}\bm{P}=\bm{0}$ and $\|(\bm{I}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime})\bm{P}\|=\zeta_{*}$. Then,

$\|(\bm{I}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime})\bm{P}\bm{P}^{\prime}\|=\|(\bm{I}-\bm{P}\bm{P}^{\prime})\hat{\bm{P}}\hat{\bm{P}}{}^{\prime}\|=\|(\bm{I}-\bm{P}\bm{P}^{\prime})\hat{\bm{P}}\|=\|(\bm{I}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime})\bm{P}\|=\zeta_{*}$** 2. 2.

$\|\bm{P}\bm{P}^{\prime}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime}\|\leq 2\|(\bm{I}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime})\bm{P}\|=2\zeta_{*}$** 3. 3.

$\|\hat{\bm{P}}{}^{\prime}\bm{Q}\|\leq\zeta_{*}$** 4. 4.

$\sqrt{1-\zeta_{*}^{2}}\leq\sigma_{i}\left((\bm{I}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime})\bm{Q}\right)\leq 1$**

Lemma 5.11**.**

[5, Lemma 3.7]** For an $n\times r$ basis matrix $\bm{P}$,

$\max_{|\mathcal{T}|\leq s}\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}\|^{2}\leq s\max_{i=1,2,\dots,r}\|\bm{I}_{i}{}^{\prime}\bm{P}\|^{2}$** 2. 2.

$\delta_{s}(\bm{I}-\bm{P}\bm{P}^{\prime})=\max_{|\mathcal{T}|\leq s}\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}\|^{2}$** 3. 3.

If $\bm{P}=[\bm{P}_{1},\bm{P}_{2}]$ then $\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}\|^{2}\leq\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}_{1}\|^{2}+\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}_{2}\|^{2}$. Also, $\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}_{1}\|^{2}\leq\|\bm{I}_{\mathcal{T}}{}^{\prime}\bm{P}\|^{2}$.

Recall that $\delta_{s}(\bm{M})$ is the $s$-restricted isometry constant [10] of a matrix $\bm{M}$, i.e., it is the smallest real number for which the following holds for all $s$-sparse vectors $\bm{z}$: $(1-\delta_{s})\|\bm{z}\|^{2}\leq\|\bm{M}\bm{z}\|^{2}\leq(1+\delta_{s})\|\bm{z}\|^{2}$.

V-C Definitions and main claim needed for Theorem 2.2 and Corollary 2.3 with ${\hat{t}}_{j}=t_{j}$

Definition 5.12**.**

We will use the following definitions in our proof.

*Let $\theta:=\theta_{j}$, $\bm{P}_{*}:=\bm{P}_{j-1}$, $\bm{P}_{\mathrm{new}}:=\bm{P}_{j,\mathrm{new}}$, $\bm{P}_{\mathrm{ch}}:=P_{j-1,{\mathrm{ch}}}$, $\bm{P}_{\mathrm{rot}}:=\bm{P}_{j,{\mathrm{rot}}}:=(\bm{P}_{{\mathrm{ch}}}\cos\theta+\bm{P}_{\mathrm{new}}\sin\theta)$, $\bm{P}:=\bm{P}_{j}$. Similarly define $\hat{\bm{P}}_{*}:=\hat{\bm{P}}_{j-1}$, $\hat{\bm{P}}:=\hat{\bm{P}}_{j}$, and *

let $\hat{\bm{P}}_{{\mathrm{rot}},k}:=\hat{\bm{P}}_{j,{\mathrm{rot}},k}$ denote the $k$-th estimate of $\bm{P}_{\mathrm{rot}}$ with $\hat{\bm{P}}_{{\mathrm{rot}},0}=[.]$. 2. 2.

*The scalars *

$\zeta_{0}^{+}:=\tilde{\varepsilon}+|\sin\theta|$*, *

$\zeta_{1}^{+}:=0.4\cdot 1.2((0.1+\tilde{\varepsilon})|\sin\theta|+\tilde{\varepsilon})+0.11\tilde{\varepsilon}$* and *

$\zeta_{k}^{+}:=0.4\cdot(1.2\zeta_{k-1}^{+})+0.11\tilde{\varepsilon}$* for $k=2,3,\dots,K$. *

We will show that these are high probability bounds on $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},k}],\bm{P}_{\mathrm{rot}})$. 3. 3.

*The events *

$\Gamma_{0}:=\{\mathrm{SE}(\hat{\bm{P}}_{*},\bm{P}_{*})\leq\tilde{\varepsilon}\}$*: clearly $\Gamma_{0}$ implies that $\mathrm{SE}(\hat{\bm{P}}_{*},\bm{P}_{\mathrm{rot}})\leq\zeta_{0}^{+}:=\tilde{\varepsilon}+|\sin\theta|$, *

$\Gamma_{k}:=\Gamma_{k-1}\cap\{\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},k}],\bm{P}_{\mathrm{rot}})\leq\zeta_{k}^{+}\}$* for $k=1,2,\dots,K$, and *

$\Gamma_{K+1}:=\Gamma_{K}\cap\{\mathrm{SE}(\hat{\bm{P}},\bm{P})\leq\tilde{\varepsilon}\}$. 4. 4.

*The time intervals: *

$\mathcal{J}_{k}:=[t_{j}+(k-1)\alpha,t_{j}+k\alpha)$* for $k=1,2,\dots,K$: the projection-SVD intervals, *

$\mathcal{J}_{K+1}:=[t_{j}+K\alpha,t_{j}+(K+1)\alpha)$*: the subspace re-estimation interval, *

$\mathcal{J}_{K+2}:=[t_{j}+K\alpha+\alpha,t_{j+1})$: the interval when the current subspace update is complete and before the next change.

We first prove the $\mathrm{SE}$ bounds of Theorem 2.2. With these, the other bounds follow easily. To obtain the $\mathrm{SE}$ bounds, we will be done if we prove the following claim.

Claim 5.13**.**

Given $\mathrm{SE}(\hat{\bm{P}}_{*},\bm{P}_{*})\leq\tilde{\varepsilon}$, w.p. at least $1-(K+1)12n^{-12}$,

$\mathrm{SE}(\hat{\bm{P}}_{*},\bm{P}_{\mathrm{rot}})\leq\zeta_{0}^{+}$, $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},1}],\bm{P}_{\mathrm{rot}})\leq\zeta_{1}^{+}$, 2. 2.

for $k>1$, $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},k}],\bm{P}_{\mathrm{rot}})\leq\zeta_{k}^{+}$, 3. 3.

and so $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},K}],\bm{P}_{\mathrm{rot}})\leq\zeta_{K}^{+}\leq\tilde{\varepsilon}$ (using definition of $K$) and $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},K}],\bm{P})\leq 2\tilde{\varepsilon}$. 4. 4.

Further, after the deletion step, $\mathrm{SE}(\hat{\bm{P}},\bm{P})\leq\tilde{\varepsilon}$.

Proving the above claim is equivalent to showing that $\Pr(\Gamma_{K+1}|\Gamma_{0})\geq 1-(K+1)12n^{-12}$.

This claim is an easy consequence of the three main lemmas and Fact 5.14 given below. Fact 5.14 provides simple upper bounds on $\zeta_{k}^{+}$ that will be used at various places. The first lemma, Lemma 5.15, shows that, assuming that the “subspace estimates so far are good enough”, the projected CS step “works” for the next $\alpha$ frames, i.e., for all $t\in\mathcal{J}_{k}$, $\hat{\mathcal{T}}_{t}=\mathcal{T}_{t}$; $\bm{e}_{t}$ is sparse and supported on $\mathcal{T}_{t}$ and satisfies (10), and $\|\bm{e}_{t}\|$ is bounded. The second lemma, Lemma 5.16, uses Lemma 5.15 and Theorem 5.7 to show that, assuming that the “subspace estimates so far are good enough”, with high probability (whp), the subspace estimate at the next projection-SVD step is even better than the previous ones. Applying Lemma 5.16 for each $k=1,2,\dots,K$ proves the first two parts of Claim 5.13. The third part follows easily from the first two and the definition of $K$. The fourth part follows using Lemma 5.17, which shows that, assuming that the $K$-th projection-SVD step produces a subspace estimate that is within $2\tilde{\varepsilon}$ of the true subspace, the subspace re-estimation step returns an estimate that is within $\tilde{\varepsilon}$ of the true subspace.

Fact 5.14**.**

Using $\tilde{\varepsilon}\leq\tilde{\varepsilon}f\leq 0.01|\sin\theta|$,

$\zeta_{0}^{+}:=\tilde{\varepsilon}+|\sin\theta|\leq 1.01|\sin\theta|$, 2. 2.

$\zeta_{1}^{+}:=0.4\cdot 1.2((0.1+\tilde{\varepsilon})|\sin\theta|+\tilde{\varepsilon})+0.11\tilde{\varepsilon}\leq 0.06|\sin\theta|$, 3. 3.

$\zeta_{k}^{+}:=0.4\cdot 1.2\zeta_{k-1}^{+}+0.11\tilde{\varepsilon}\leq 0.5^{k-1}\zeta_{1}^{+}+\frac{0.11}{1-0.5}\tilde{\varepsilon}\leq 0.5^{k-1}(0.06|\sin\theta|)+0.11\tilde{\varepsilon}\leq 0.03|\sin\theta|$

This claim essentially implies Theorem 2.2 and Corollary 2.3 with ${\hat{t}}_{j}=t_{j}$. We prove these without this assumption in Appendix A.

V-D The three main lemmas needed to prove the main claim and their proofs

Lemma 5.15** (Projected CS).**

Recall from Sec. V that $s$ is an upper bound on $|\mathcal{T}_{t}|$. Under assumptions of Theorem 2.2 or or Corollary 2.3, the following hold for $k=1,2,\dots,K+2$. Let $\bm{\Psi}_{1}:=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}$, $\bm{\Psi}_{k}:=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}-\hat{\bm{P}}_{{\mathrm{rot}},k-1}\hat{\bm{P}}_{{\mathrm{rot}},k-1}{}^{\prime}$ for $k=2,3,\dots,K+1$, and $\bm{\Psi}_{K+2}:=\bm{I}-\hat{\bm{P}}\hat{\bm{P}}^{\prime}$. From Algorithm 1,

[TABLE]

Assume that $\Gamma_{k-1}$ holds. Then,

$\max_{|\mathcal{T}|\leq 2s}\|\bm{I}_{\mathcal{T}}{}^{\prime}\hat{\bm{P}}_{*}\|\leq 0.3+\tilde{\varepsilon}\leq 0.31$. 2. 2.

$\max_{|\mathcal{T}|\leq 2s}\|\bm{I}_{\mathcal{T}}{}^{\prime}\hat{\bm{P}}_{{\mathrm{rot}},k-1}\|\leq 0.1+\tilde{\varepsilon}+\frac{\zeta_{k-1}^{+}+\tilde{\varepsilon}}{|\sin\theta|}\leq 0.1+0.01+0.04<0.15$. 3. 3.

$\delta_{2s}(\bm{\Psi}_{1})\leq 0.31^{2}<0.12$, $\delta_{2s}(\bm{\Psi}_{k})\leq 0.31^{2}+0.15^{2}<0.12$ for $k=2,3,\dots,K+2$, 4. 4.

for all $t\in\mathcal{J}_{k}$, $\|(\bm{\Psi}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}_{\mathcal{T}_{t}})^{-1}\|\leq 1.2$ 5. 5.

for all $t\in\mathcal{J}_{k}$, $\hat{\mathcal{T}}_{t}=\mathcal{T}_{t}$ 6. 6.

for all $t\in\mathcal{J}_{k}$, $\bm{e}_{t}:=\bm{\hat{x}}_{t}-\bm{x}_{t}=\bm{\ell}_{t}-\hat{\bm{\ell}}_{t}+\bm{v}_{t}$ satisfies (10) with $\bm{\Psi}=\bm{\Psi}_{k}$ for $t\in\mathcal{J}_{k}$ 7. 7.

*for $t\in\mathcal{J}_{1}$, $\|\bm{e}_{t}\|\leq 2.4\sqrt{\eta}(\tilde{\varepsilon}\sqrt{r\lambda^{+}}+0.11\Delta\sqrt{\lambda_{\mathrm{ch}}})$; *

*for $t\in\mathcal{J}_{k}$, $\|\bm{e}_{t}\|\leq 2.4\sqrt{\eta}(\tilde{\varepsilon}\sqrt{r\lambda^{+}}+\zeta_{k-1}^{+}\sqrt{\lambda_{\mathrm{ch}}})$ for $k=2,3,\dots K$; *

*for $t\in\mathcal{J}_{K+1}$, $\|\bm{e}_{t}\|\leq 4.8\sqrt{\eta}\tilde{\varepsilon}\sqrt{r\lambda^{+}}$. *

for $t\in\mathcal{J}_{K+2}$, $\|\bm{e}_{t}\|\leq 2.4\sqrt{\eta}\tilde{\varepsilon}\sqrt{r\lambda^{+}}$.

Proof.

Since the noise bound of Theorem 2.2 is much smaller or equal to those assumed by Corollary 2.3, if we can prove the latter, we would have also proved the former. Using the first claim of Lemma 5.11, the max-outlier-frac-col bound of Theorem 2.2 and the incoherence assumption (6) imply that, for any set $\mathcal{T}$ with $|\mathcal{T}|\leq 2s$,

[TABLE]

(In order to simplify our assumptions, we have simplified the incoherence/denseness assumption from what it was in the original version of this work; as a result, even the first term above is bounded by 0.1 (not 0.3 as before). To not have to change the rest of the proof given below, we still use the 0.3 bound in the writing below). Using (13), for any set $\mathcal{T}$ with $|\mathcal{T}|\leq 2s$,

[TABLE]

The second row used (13) and the following: since $\hat{\bm{P}}_{*}$ and $\bm{P}_{*}$ have the same dimension, $\mathrm{SE}(\bm{P}_{*},\hat{\bm{P}}_{*})=\mathrm{SE}(\hat{\bm{P}}_{*},\bm{P}_{*})$ (follows by Lemma 5.10, item 1). The third row follows using the definition of event $\Gamma_{k-1}$. Proceeding similarly for $\hat{\bm{P}}_{{\mathrm{rot}},k-1}$ and $\bm{P}_{\mathrm{new}}$ (both have the same dimension),

[TABLE]

The second row used item 1 of Lemma 5.10 and (13). The third row used triangle inequality. The last row follows using $\bm{P}_{\mathrm{new}}=\frac{\bm{P}_{\mathrm{rot}}-\bm{P}_{{\mathrm{ch}}}\cos\theta}{\sin\theta}$ and the definition of event $\Gamma_{k-1}$. Using this, triangle inequality and $|\cos\theta|\leq 1$, we bound the first term. Using item 3 of Lemma 5.10, we bound the second term. The final bounds use Fact 5.14.

To get the above bound, we use $\bm{P}_{\mathrm{new}}$ (and not $\bm{P}_{\mathrm{rot}}$) because $\|\hat{\bm{P}}_{*}{}^{\prime}\bm{P}_{\mathrm{new}}\|\leq\tilde{\varepsilon}$ since $\bm{P}_{*}$ is orthogonal to $\bm{P}_{\mathrm{new}}$. But we do not have a small upper bound on $\|\hat{\bm{P}}_{*}{}^{\prime}\bm{P}_{\mathrm{rot}}\|$.

The third claim follows using the first two claims and Lemma 5.11. The fourth claim follows from the third claim as follows:

[TABLE]

The last three claims follow the approach of the proof of [5, Lemma 6.4]. There are minor differences because we set $\xi$ a little differently now and because we assume $\bm{v}_{t}\neq 0$. We provide the proof in Appendix C. ∎

Lemma 5.16** (Projection-SVD).**

Under the assumptions of Theorem 2.2 or Corollary 2.3, the following holds for $k=1,2,\dots K$. Conditioned on $\Gamma_{k-1}$, w.p. at least $1-12n^{-12}$, $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},k}],\bm{P}_{\mathrm{rot}})\leq\zeta_{k}^{+}$, i.e., $\Gamma_{k}$ holds.

Proof.

Since the noise bound of Theorem 2.2 is much smaller or equal to those assumed by Corollary 2.3, if we can prove the latter, we would have also proved the former.

Assume that $\Gamma_{k-1}$ holds. The proof first uses Lemma 5.15 to get an expression for $\bm{e}_{t}=\bm{\ell}_{t}-\hat{\bm{\ell}}_{t}+\bm{v}_{t}$ and then applies Theorem 5.7 with the modification given in Remark 5.8. Using Lemma 5.15, for all $t\in\mathcal{J}_{k}$,

[TABLE]

where $\bm{\Psi}=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}-\hat{\bm{P}}_{{\mathrm{rot}},k-1}\hat{\bm{P}}_{{\mathrm{rot}},k-1}{}^{\prime}$ with $\hat{\bm{P}}_{{\mathrm{rot}},0}=[.]$.

In the $k$-th projection-SVD step, we use these $\hat{\bm{\ell}}_{t}$’s and $\hat{\bm{P}}_{*}$ to get a new estimate of $\bm{P}_{\mathrm{rot}}$ using projection-SVD. To bound $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},k}],\bm{P}_{\mathrm{rot}})$, we apply Theorem 5.7 (Remark 5.8)111111We use Remark 5.8 with $\mathcal{E}_{*}\equiv\Gamma_{0}$, $Z\equiv\{\hat{\bm{\ell}}_{1},\hat{\bm{\ell}}_{2},\dots,\hat{\bm{\ell}}_{t_{j}+(k-1)\alpha-1}\}$, and $\tilde{\mathcal{E}}(Z)=\Gamma_{k-1}\setminus\Gamma_{0}$. with $\mathcal{E}_{0}\equiv\Gamma_{k-1}$, $\bm{y}_{t}\equiv\hat{\bm{\ell}}_{t}$, $\bm{w}_{t}\equiv-\bm{e}_{l,t}$, $\bm{z}_{t}\equiv-\bm{e}_{v,t}+\bm{v}_{t}$, $\alpha\geq\alpha_{0}\equiv\alpha_{*}$, and $\mathcal{J}^{\alpha}\equiv\mathcal{J}_{k}$. We can let $\bm{M}_{2,t}=-\bm{I}_{\mathcal{T}_{t}}$ which implies $b_{0}\equiv\text{\small{max-outlier-frac-row}}^{\alpha}$ and $\bm{M}_{1,t}=(\bm{\Psi}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}_{\mathcal{T}_{t}})^{-1}\bm{I}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}$. Using the $\text{\small{max-outlier-frac-row}}^{\alpha}$ bound of Theorem 2.2 and Lemma 5.9, the main assumption needed by Theorem 5.7, (11), holds. With $\bm{P}=\bm{P}_{j}$ satisfying (5), and $\alpha_{*}$ defined in (9), all the key assumptions of Theorem 5.7 hold. The simpler expression of $\alpha_{*}$ suffices because we treat $\eta$ as a numerical constant and so $f^{2}(r\log n)>f^{2}(r+\log n)$ for large $n,r$.

We now just need to compute $q_{0}$ and $q_{\mathrm{rot}}$ for each $k$, ensure that they satisfy (12), and apply the result. The computation for $k=1$ is different from the rest. When $k=1$, $\bm{\Psi}=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}$. Thus, using item 4 of Lemma 5.15 and the definition of event $\Gamma_{k-1}$, $\|\bm{M}_{1,t}\bm{P}_{*}\|\leq 1.2\tilde{\varepsilon}=q_{0}$, $q_{0}<2\tilde{\varepsilon}$, and

[TABLE]

The third row follows using (13) and $\|\hat{\bm{P}}_{*}{}^{\prime}\bm{P}_{\mathrm{new}}\|\leq\tilde{\varepsilon}$ (folows by item 3 of Lemma 5.10). Using $\tilde{\varepsilon}\leq\tilde{\varepsilon}f\leq 0.01|\sin\theta|$, clearly $q_{\mathrm{rot}}<0.2|\sin\theta|$. Finally, in this interval, the bound on $b_{z}$ is satisfied since $b_{z}=b_{v,t}$ and the expression for $b_{v,t}$ in $\mathcal{J}_{1}$ is equal to the required upper bound on $b_{z}$. Applying Theorem 5.7, $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},1}],\bm{P}_{\mathrm{rot}})\leq 0.4q_{\mathrm{rot}}+0.11\tilde{\varepsilon}=0.4\cdot 1.2((0.1+\tilde{\varepsilon})|\sin\theta|+\tilde{\varepsilon})+0.11\tilde{\varepsilon}=\zeta_{1}^{+}$.

Consider $k>1$. Now $\bm{\Psi}=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}-\hat{\bm{P}}_{{\mathrm{rot}},k-1}\hat{\bm{P}}_{{\mathrm{rot}},k-1}{}^{\prime}$. With this, we still have $\|\bm{M}_{1,t}\bm{P}_{*}\|\leq 1.2\tilde{\varepsilon}=q_{0}$ and $q_{0}<2\tilde{\varepsilon}$. But, to bound $\|\bm{M}_{1,t}\bm{P}_{\mathrm{rot}}\|$ we cannot use the approach that worked for $k=1$. The reason is that $\|[\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},k-1}]^{\prime}\bm{P}_{\mathrm{new}}\|$ is not small. However, instead, we can now use the fact that $[\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},k-1}]$ is a good estimate of $\bm{P}_{\mathrm{rot}}$, with $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},k-1}],\bm{P}_{\mathrm{rot}})\leq\zeta_{k-1}^{+}$ (from definition of event $\Gamma_{k-1}$). Thus,

[TABLE]

By Fact 5.14, $q_{\mathrm{rot}}<0.2|\sin\theta|$. Even in this interval, the required bound on $b_{z}$ holds. Thus, applying Theorem 5.7, $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},k}],\bm{P}_{\mathrm{rot}})\leq 0.4q_{\mathrm{rot}}+0.11\tilde{\varepsilon}=0.4\cdot 1.2\zeta_{k-1}^{+}+0.11\tilde{\varepsilon}=\zeta_{k}^{+}$. ∎

Lemma 5.17** (Simple SVD based subspace re-estimation).**

Under the assumptions of Theorem 2.2 or Corollary 2.3, the following holds. Conditioned on $\Gamma_{K}$, w.p. at least $1-12n^{-12}$, $\mathrm{SE}(\hat{\bm{P}},\bm{P})\leq\tilde{\varepsilon}$, i.e., $\Gamma_{K+1}$ holds.

Proof.

Assume that $\Gamma_{K}$ holds. Using Lemma 5.15, for all $t\in\mathcal{J}_{K+1}$,

[TABLE]

where $\bm{\Psi}=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}-\hat{\bm{P}}_{{\mathrm{rot}},K}\hat{\bm{P}}_{{\mathrm{rot}},K}{}^{\prime}$. Re-estimating the entire subspace using simple SVD applied to these $\hat{\bm{\ell}}_{t}$’s is an instance of correlated-PCA with $\bm{y}_{t}\equiv\hat{\bm{\ell}}_{t}$, $\bm{w}_{t}\equiv-\bm{e}_{l,t}$ and $\bm{z}_{t}\equiv-\bm{e}_{v,t}+\bm{v}_{t}$. We can apply the following result for correlated-PCA [23, Theorem 2.13] to bound $\mathrm{SE}(\hat{\bm{P}},\bm{P})$. Recall $\hat{\bm{P}}$ contains the top $r$ eigenvectors of $\sum_{t\in\mathcal{J}_{K+1}}\hat{\bm{\ell}}_{t}\hat{\bm{\ell}}_{t}^{\prime}$. The following is a simplified version of [23, Theorem 2.13]. It follows by upper bounding $\lambda_{z,\bm{P},\bm{P}_{\perp}}$ and $\lambda_{z,rest}^{+}$ by $\lambda_{z}^{+}$ and lower bound $\lambda_{z,\bm{P}}^{-}$ by zero in [23, Theorem 2.13].

Theorem 5.18**.**

For $t\in\mathcal{J}^{\alpha}$, we are given data vectors $\bm{y}_{t}:=\bm{\ell}_{t}+\bm{w}_{t}+\bm{z}_{t}$ where $\bm{w}_{t}=\bm{M}_{t}\bm{\ell}_{t}$, $\bm{\ell}_{t}=\bm{P}\bm{a}_{t}$ and $\bm{z}_{t}$ is small unstructured noise. Let $\hat{\bm{P}}$ be the matrix of top $r$ eigenvectors of $\frac{1}{\alpha}\sum_{t\in\mathcal{J}^{\alpha}}\bm{y}_{t}\bm{y}_{t}{}^{\prime}$. Assume that $\bm{M}_{t}$ can be decomposed as $\bm{M}_{t}=\bm{M}_{2,t}\bm{M}_{1,t}$ so that $\|\bm{M}_{2,t}\|\leq 1$ but $\left\|\frac{1}{\alpha}\sum_{t}\bm{M}_{2,t}\bm{M}_{2,t}{}^{\prime}\right\|\leq b$ for a $b<1$. Let $q$ be an upper bound on $\max_{t\in\mathcal{J}^{\alpha}}\|\bm{M}_{1,t}\bm{P}\|$. We assume that $\|\bm{z}_{t}\|\leq b_{z}$ and define $\|\mathbb{E}{[\bm{z}_{t}\bm{z}_{t}{}^{\prime}]}\|\leq\lambda_{z}^{+}:=b_{z}^{2}/r$. For an $\varepsilon_{\mathrm{SE}}>0$, define

[TABLE]

If $\alpha\geq\alpha_{0}$, and $3\sqrt{b}qf+\frac{\lambda_{z}^{+}}{\lambda^{-}}\leq 0.46\varepsilon_{\mathrm{SE}},$ then, w.p. at least $1-12n^{-12}$, $\mathrm{SE}(\hat{\bm{P}},\bm{P})\leq\varepsilon_{\mathrm{SE}}.$ .

Apply the above result with $\bm{y}_{t}\equiv\hat{\bm{\ell}}_{t}$, $\bm{w}_{t}\equiv-\bm{e}_{l,t}$, $\bm{z}_{t}\equiv-\bm{e}_{v,t}+\bm{v}_{t}$, $\alpha\geq\alpha_{*}$, and $\mathcal{J}^{\alpha}\equiv\mathcal{J}_{K+1}$. From the expression for $\bm{e}_{t}$, we can let $\bm{M}_{2,t}\equiv-\bm{I}_{\mathcal{T}_{t}}$, $\bm{M}_{1,t}\equiv(\bm{\Psi}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}_{\mathcal{T}_{t}})^{-1}\bm{I}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}$. Next we compute $q$. Since $\Gamma_{K}$ holds, $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},K}],\bm{P})\leq\tilde{\varepsilon}+\zeta_{K}^{+}\leq 2\tilde{\varepsilon}$. Thus, $\|\bm{M}_{1,t}\bm{P}\|\leq 1.2\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{{\mathrm{rot}},K}],\bm{P})\leq 1.2\cdot 2\tilde{\varepsilon}=q$. The final desired error is $\varepsilon_{\mathrm{SE}}=\tilde{\varepsilon}$. Using Lemma 5.9 and the $\text{\small{max-outlier-frac-row}}^{\alpha}$ bound from Theorem 2.2, the bound on the time-average of $\bm{M}_{2,t}\bm{M}_{2,t}{}^{\prime}$ holds with $b\equiv\rho_{\mathrm{row}}=\frac{0.01}{f^{2}}<\frac{0.5^{2}}{(3\cdot 2.4f)^{2}}$. Also, in this interval $b_{z}^{2}=Cr\lambda^{+}$ and so $\lambda_{z}=C\lambda^{+}$ and thus the second term in the $\alpha_{0}$ expression equals the first term. The third term can of course be ignored in the large $n,r$ regime. Applying the above result with $\varepsilon_{\mathrm{SE}}=\tilde{\varepsilon}$, and $q=2.4\tilde{\varepsilon}$, we conclude the following: for $\alpha\geq\alpha_{*}$, w.p. at least $1-12n^{-12}$, $\mathrm{SE}(\hat{\bm{P}},\bm{P})\leq\tilde{\varepsilon}$. The simpler expression of $\alpha_{*}$ suffices because $\eta$ is treated as a numerical constant and so $f^{2}(r\log n)>f^{2}(r+\log n)$ for large $n,r$. Also under the assumption of Corollary 2.3, the second term is dominated by the first term. ∎

Proof of Claim 5.13.

Lemma 5.16 tells us that $\Pr(\Gamma_{k}|\Gamma_{k-1})\geq 1-12n^{-12}$. Lemma 5.17 tells us that $\Pr(\Gamma_{K+1}|\Gamma_{K})\geq 1-12n^{-12}$. Thus, $\Pr(\Gamma_{K+1}|\Gamma_{0})=\Pr(\Gamma_{K+1},\Gamma_{K},\dots\Gamma_{1}|\Gamma_{0})=\Pr(\Gamma_{1}|\Gamma_{0})\Pr(\Gamma_{2}|\Gamma_{1})\dots\Pr(\Gamma_{K+1}|\Gamma_{K})\geq(1-12n^{-12})^{K}(1-12n^{-12}).$ since $\Gamma_{K+1}\subseteq\Gamma_{K}\subseteq\Gamma_{K-1}\dots\subseteq\Gamma_{0}$. The result follows since $(1-12n^{-12})^{K}(1-12n^{-12})\geq 1-(K+1)12n^{-12}$. ∎

Proof of Theorem 2.2 with ${\hat{t}}_{j}=t_{j}$.

Define the events $\Gamma_{1,0}:=\{\mathrm{SE}(\hat{\bm{P}}_{0},\bm{P}_{0})\leq\tilde{\varepsilon}\}$, $\Gamma_{j,k}:=\Gamma_{j,k-1}\cap\{\mathrm{SE}([\hat{\bm{P}}_{j-1},\hat{\bm{P}}_{j,{\mathrm{rot}},k}])\leq\zeta_{k}^{+}\}$, for $k=1,2,\dots,K$, $\Gamma_{j,K+1}:=\Gamma_{j,K}\cap\{\mathrm{SE}(\hat{\bm{P}}_{j},\bm{P}_{j})\leq\tilde{\varepsilon}\}$ and $\Gamma_{j+1,0}:=\Gamma_{j,K+1}$. We can state and prove Lemmas 5.16 and 5.17 with $\Gamma_{k}$ replaced by $\Gamma_{j,k}$. Then Claim 5.13 implies that $\Pr(\Gamma_{j,K+1}|\Gamma_{j,0})\geq 1-12n^{-12}$. Using $\Gamma_{J,K+1}\subseteq\Gamma_{J-1,K+1}\dots\subseteq\Gamma_{1,K+1}\subseteq\Gamma_{1,0}$ and $\Gamma_{j+1,0}:=\Gamma_{j,K+1}$, $\Pr(\Gamma_{J,K+1}|\Gamma_{1,0})=\Pr(\Gamma_{J,K+1},\Gamma_{J-1,K+1},\dots\Gamma_{1,K+1}|\Gamma_{1,0})=\Pr(\Gamma_{1,K+1}|\Gamma_{1,0})\Pr(\Gamma_{2,K+1}|\Gamma_{2,0})\dots\Pr(\Gamma_{J,K+1}|\Gamma_{J,0})\geq(1-(K+1)12n^{-12})^{J}\geq 1-J(K+1)12n^{-12}\geq 1-dn^{-12}$.

Event $\Gamma_{J,K+1}$ implies that $\Gamma_{j,k}$ holds for all $j$ and for all $k$. Thus, all the $\mathrm{SE}$ bounds given in Theorem 2.2 hold. Using Lemma 5.15, $\hat{\mathcal{T}}_{t}=\mathcal{T}_{t}$ for all the time intervals of interest, and the bounds on $\|\bm{e}_{t}\|$ hold. ∎

VI Empirical Evaluation

In this section we illustrate the superiority of s-ReProCS over existing state of the art methods on synthetic and real data. In particular, we consider the task of background subtraction. All time comparisons are performed on a Desktop Computer with Intel${}^{\textsuperscript{\textregistered}}$ Xeon E$3$-$1240$ $8$-core CPU @ $3.50$GHz and $32$GB RAM. And all experiments with synthetic data are averaged over $100$ independent trials. All codes are available at https://github.com/praneethmurthy/ReProCS.

Similar experiments have been shown in the earlier ReProCS works (original-ReProCS) [5, 9, 12, 7]. The purpose of this section is to illustrate that, even though s-ReProCS is much simpler, is provably faster and memory efficient, and provably works under much simpler assumptions, its experimental performance is still similar to that of original-ReProCS. It outperforms existing works for the same classes of videos and simulated data for which original-reprocs outperforms them.

VI-A Synthetic Data

Our first simulation experiment is done to illustrate the advantage of s-ReProCS over existing batch and online RPCA techniques. As explained earlier, because s-ReProCS exploits dynamics (slow subspace change), it is provably able to tolerate a much larger fraction of outliers per row than all the existing techniques without needing uniformly randomly generated support sets. When the number of subspace changes, $J$, is large, it also tolerates a significantly larger fraction of outliers per column. The latter is hard to demonstrate via simulations (making $J$ large will require a very long sequence). Thus we demonstrate only the former. Our second experiment shows results with using an i.i.d. Bernoulli model on support change (which is the model assumed in the other works).

One practical instance where outlier fractions per row can be larger than those per column is in the case of video moving objects that are either occasionally static or slow moving [12, 7]. The outlier support model for our first and second experiments is inspired by this example and the model used in [12, 7]. It models a 1D video consisting of a person/object of length $s$ pacing up and down in a room with frequent stops. The object is static for $\beta$ frames at a time and then moves down. It keeps moving down for a period of $\tau$ frames, after which it turns back up and does the same thing in the other direction. We let $\beta=\lceil c_{0}\tau\rceil$ for a $c_{0}<1$. With this model, for any interval of the form $[(k_{1}-1)\tau+1,k_{2}\tau]$ for $k_{1},k_{2}$ integers, the outlier fraction per row is bounded by $c_{0}$. For any general interval of length $\alpha\geq\tau$, this $\text{\small{max-outlier-frac-row}}^{\alpha}$ is still bounded by $2c_{0}$ while max-outlier-frac-col is bounded by $s/n$.

Model 6.19**.**

Let $\beta=\lceil c_{0}\tau\rceil$. Assume that the $\mathcal{T}_{t}$ satisfies the following. For the first $\tau$ frames (downward motion),

[TABLE]

for the next $\tau$ frames (upward motion), $\mathcal{T}_{t}=$

[TABLE]

Starting at $t=2\tau+1$, the above pattern is repeated every $2\tau$ frames until the end, $t=d$.

This model is motivated by the model assumed for the guarantees in older works [12, 7]. The above model is one practically motivated way to simulate data that is not not generated uniformly at random (or as i.i.d. Bernoulli, which is approximately the same as the uniform model for large $n$). It also provides a way to generate data with a different bounds on outlier fractions per row and per column. The maximum outlier fraction per column is $s/n$. For any time interval of length $\alpha\geq\tau$, the outlier fraction per row is bounded by $2c_{0}$. Thus, for Theorem 2.2, with this model, $\rho_{\mathrm{row}}=2c_{0}/f^{2}$. By picking $2c_{0}$ larger than $s/n$ we can ensure larger outlier fractions per row than per column.

We compare Algorithm 4 and its offline counterpart with three of the batch methods with provably guarantees discussed in Sec. II-C - PCP [2], AltProj [6] and RPCA-GD [8] - and with two recently proposed online algorithms known to have good experimental performance and for which code was available - ORPCA [14] and GRASTA [15]. The code for all these techniques are cloned from the Low-Rank and Sparse library (https://github.com/andrewssobral/lrslibrary).

For generating data we used $d=8000$, $t_{\mathrm{train}}=500$, $J=2$, $r=5$, $f=16$ with $t_{1}=1000$, $\theta_{1}=30^{\circ}$, $t_{2}=4300$, $\theta_{2}=1.01\theta_{1}$ and varying $n$. The $\bm{a}_{t}$’s are zero mean i.i.d uniform random variables generated exactly as described before and so is $\bm{P}_{(t)}$. We generated a basis matrix $\bm{Q}$ by ortho-normalizing the first $r+2$ columns of a $n\times n$ i.i.d. standard Gaussian matrix. For $t\in[1,\ t_{1})$, we set $\bm{P}_{(t)}=\bm{P}_{0}$ with $\bm{P}_{0}$ being the first $r$ columns of $\bm{Q}$. We let $\bm{P}_{1,\mathrm{new}}$ be $(r+1)$-th column of $\bm{Q}$, and rotated it in using (5) with $\bm{U}_{1}=\bm{I}$ and with angle $\theta_{1}$ to get $\bm{P}_{1}$. We set $\bm{P}_{(t)}=\bm{P}_{1}$ for $t\in[t_{1},t_{2})$. We set $\bm{P}_{2,\mathrm{new}}$ to be the last column of $\bm{Q}$, $\bm{U}_{2}=\bm{I}$, and rotate using angle $\theta_{2}$ just as done in the first subspace change and finally for $t\in[t_{2},d]$ we set $\bm{P}_{(t)}=\bm{P}_{2}$. At all times $t$, we let $\bm{\ell}_{t}=\bm{P}_{(t)}\bm{a}_{t}$ with $\bm{a}_{t}$ being zero mean, i.i.d uniform random variables such that $(\bm{a}_{t})_{i}\sim unif\left(-\sqrt{f},\ \sqrt{f}\right)$ for $i=1,\ \cdots,\ r-2$ and $(\bm{a}_{t})_{[r-1,r]}\sim unif(-1,1)$. With this the condition number is $f$, the covariance matrix, $\Lambda=\mathrm{diag}(f,f,\cdots,f,1,1)/3$, $\lambda^{+}=f/3$, $\lambda_{\mathrm{ch}}=\lambda^{-}=1/3$, and $\eta=3$. We generate $\mathcal{T}_{t}$ using Model 6.19 as follows. For $t\in[t_{\mathrm{train}},d]$, we used $s=0.1n$, $c_{0}=0.2$ and $\tau=100$. Thus $\rho_{\mathrm{row}}=0.4/f^{2}$. For $t\in[1,t_{\mathrm{train}}]$, we used $s=0.05n$ and $c_{0}=0.02$. This was done to ensure that AltProj (or any other batch technique works well for this period and provides a good initialization). The magnitudes of the nonzero entries of $\bm{x}_{t}$ (outliers) were generated s i.i.d uniform r.v.’s between $x_{\min}=10$ and $x_{\max}=25$.

We implemented Algorithm 4 for s-ReProCS (with initialization using AltProj) with $\alpha=Cf^{2}r\log n=500$, $K=\lceil-0.8\log(0.9\tilde{\varepsilon})\rceil=5$, $\omega_{supp}=x_{\min}/2$, and $\omega_{evals}=0.0025\lambda^{-}$. We initialized using AltProj applied to $\bm{Y}_{[1,t_{\mathrm{train}}]}$. For the batch methods used in the comparisons – PCP, AltProj and RPCA-GD, we implement the algorithms on $\bm{Y}_{[1,t]}$ every $t=t_{\mathrm{train}}+k\alpha-1$ frames. Further, we set the regularization parameter for PCP $\lambda=1/\sqrt{n}$ in accordance with [2]. The other known parameters, $r$ for Alt-Proj, outlier-fraction for RPCA-GD, are set using the true values. For online methods we implement the algorithms without modifications. The regularization parameter for ORPCA was set as with $\lambda_{1}=1/\sqrt{n}$ and $\lambda_{2}=1/\sqrt{d}$ according to [14]. We plot the subspace error and the normalized error of $\bm{\ell}_{t}$ over time in Fig. 2(a) and 2(b) for $n=5000$. We display the time-averaged error for other values of $n$ in Table V. This table also contains the time comparisons.

As can be seen, s-ReProCS outperforms all the other methods and offline s-ReProCS significantly outperforms all the other methods for this experiment. The reason is that the outlier fraction per row are quite large, but s-ReProCS exploits slow subspace change. In principle, even GRASTA exploits slow subspace change, however, it uses approximate methods for computing the SVD and does not use projection-SVD and hence it fails. s-ReProCS and offline s-ReProCS are faster than all the batch methods especially for large $n$. In fact when $n=10000$, the batch methods are out of memory and cannot work, while s-ReProCS still can. But s-ReProCS is slower than GRASTA and ORPCA.

*Comparison with other algorithms - random outlier support using the i.i.d. Bernoulli model. * We generated data exactly as described above with the following change: $\mathcal{T}_{t}$ was now generated as i.i.d. Bernoulli with probability of any index $i$ being in $\cup_{t\in[1,n]}\mathcal{T}_{t}$ being $\rho_{s}=0.02$ for the first $t_{\mathrm{train}}$ frames and $\rho_{s}=0.2$ for the subsequent data. Notice that under the Bernoulli model, $\rho_{\mathrm{row}}=\rho_{\mathrm{col}}=\rho_{s}$. We used $n=500$. We show the results in Fig. 2(c) and 2(d). For this experiment, the batch methods PCP and AltProj have good performance, that is better than s-ReProCS at most time instants. Offline s-ReProCS still outperforms all the other methods.

VI-B Real Data: Background Subtraction

In this section we provide simulation results for on real videos on three benchmark datasets. For all the sequences, to implement s-ReProCS, we obtained an estimate using the AltProj algorithm. For the initialization we set $r=40$ and the other parameters for the proposed algorithm, were set as follows. We used $\alpha=60$, $K=3$, $\xi_{t}=\|\bm{\Psi}\hat{\bm{\ell}}_{t-1}\|_{2}$ and $\omega_{evals}=0.0011\lambda^{-}$. We found that these parameters work for most videos that we verified our algorithm on. For a more detailed empirical evaluation on real world data-sets, please see [21]. The other state-of-the-art algorithms were implemented using the default setting. In algorithms where we are required to provide an esimate of the rank, we used $r=40$ consistently. Additionally, for RPCA-GD we set the corruption fraction, $\alpha=0.2$ as described in the paper. We must also mention that the performance of s-ReProCS w.r.t. original-ReProCS is very similar, but is provably fast, and needs fewer assumptions. Again, a more detailed comparison is presented in [21].

Meeting Room (MR) dataset: The meeting room sequence is a set of $1964$ images of resolution $64\times 80$. The first $1755$ frames consists of outlier-free data and so we only consider the last $1209$ frames. Here and below, we use the first $400$ “noisy” frames as the training data and the algorithm parameters are set as mentioned before. This is a challenging video sequence because the color of the person and the color of the curtain are hard to distinguish. s-ReProCS algorithm is able to perform the separation at around $43$ frames per second. The recovered background images are shown in the first two rows of Fig. 3.

Switch Light (SL) dataset: This dataset contains $2100$ images of resolution $120\times 160$. The first $770$ frames are outlier free. This is a challenging sequence because there are drastic changes in the subspace as indicated in the last two rows of Fig. 3. This causes all the batch techniques to fail. For this sequence, s-ReProCS achieves a “test” processing rate of $16$ frames-per-second. The recovered background images are shown in the middle two rows of Fig. 3.

Lobby (LB) dataset: This dataset contains $1555$ images of resolution $128\times 160$. The first $341$ frames are outlier free. This is a challenging sequence, as the background changes often due to illumination changes, and there are multiple objects in the foreground to detect and subtract. For this sequence, s-ReProCS achieves a “test” processing rate of $12$ frames-per-second. The images are shown in the last two rows of Fig. 3.

VII Conclusions and Future Work

We obtained the first complete guarantee for any online, streaming or dynamic RPCA algorithm that holds under weakened versions of standard RPCA assumptions, slow subspace change, and outlier magnitudes are either large or very small. Our guarantee implies that, by exploiting these extra assumptions, one can significantly weaken the required bound on outlier fractions per row. This has many important practical implications especially for video analytics. We analyzed a simple algorithm based on the Recursive Projected Compressive Sensing (ReProCS) framework introduced in [5]. The algorithm itself is simpler than other previously studied ReProCS-based methods, it is provably faster, and has near-optimal memory complexity. Moreover, our guarantee removes all the strong assumptions made by the previous two guarantees for ReProCS-based methods.

As described earlier, our current result still has limitations, some of which can be removed with a little more work. For example, it assumes a very simple model on subspace change in which only one direction can change at any given change time. Of course the changing direction could be different at different change times, and hence over a long period, the entire subspace could change. In follow-up work [11], we have studied another algorithm that removes this limitation. Another issue that we would like to study is whether the lower bound on outlier magnitudes can be relaxed further if we use the stronger assumption on outlier fractions per row (assume they are of order $1/r$). It may be possible to do this by borrowing the AltProj [6] proof idea.

A question of practical and theoretical interest is to develop a streaming version of simple-ReProCS for dynamic RPCA. A preprint that studies a streaming algorithm for standard RPCA but only for the restrictive $r_{\scriptscriptstyle{L}}=r=1$ setting is [34]. By streaming we mean that the algorithm makes only one pass through the data and needs storage of order exactly $nr$. Simple-ReProCS needs only a little more storage than this, however, it makes multiple passes through the data in the SVD steps. Algorithmically, streaming ReProCS is easy to develop: one can replace the projection SVD and SVD steps in the subspace update by their streaming versions, e.g., block stochastic power method. However, in order to prove that this still works (with maybe an extra factor of $\log n$ in the delay), one would need to analyze the block stochastic power method for the problems of PCA in data-dependent noise, and for its extension that assumes availability of partial subspace knowledge. Finally, as explained in [12], any guarantee for dynamic RPCA also provides a guarantee for dynamic Matrix Completion (MC) as an almost direct corollary. The reason is that MC can be interpreted as RPCA with outlier supports $\mathcal{T}_{t}$ being known.

Appendix A Proof of Theorem 2.2 or Corollary 2.3 without assuming $t_{j}$ known

The key results needed for this and later proofs – Cauchy-Schwarz for sums of matrices, matrix Bernstein, and Vershynin’s sub-Gaussian result – are summarized in the last appendix, Appendix E. For a summary of notation used for this and later proofs, please see Table IV.

Here we prove Theorem 2.2 in the general case. The main idea is explained in Sec. IV-B. Define

[TABLE]

Thus, ${\hat{t}}_{j-1,fin}$ is the time at which the $(j-1)$-th subspace update is complete; whp, this occurs before $t_{j}$. Under this assumption, $t_{j,*}$ is such that $t_{j}$ lies in the interval $[t_{j,*}-\alpha+1,t_{j,*}]$. Recall from the algorithm that we increment $j$ to $j+1$ at $t={\hat{t}}_{j}+K\alpha+\alpha:={\hat{t}}_{j,fin}$. Thus, for $t\in[t_{j},{\hat{t}}_{j,fin})$, $\bm{\Phi}=\bm{I}-\hat{\bm{P}}_{\ast}\hat{\bm{P}}_{\ast}{}^{\prime}$, while for $t\in[{\hat{t}}_{j,fin},t_{j+1})$, $\bm{\Phi}=\bm{I}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime}$.

Definition A.20**.**

Define the events

$\mathrm{Det0}:=\{{\hat{t}}_{j}=t_{j,*}\}=\{\lambda_{\max}(\frac{1}{\alpha}\sum_{t=t_{j,*}-\alpha+1}^{t_{j,*}}(\bm{I}-\hat{\bm{P}}_{\ast}\hat{\bm{P}}_{\ast}{}^{\prime})\hat{\bm{\ell}}_{t}\hat{\bm{\ell}}_{t}^{\prime}(\bm{I}-\hat{\bm{P}}_{\ast}\hat{\bm{P}}_{\ast}{}^{\prime}))>\omega_{evals}\}$* and *

$\mathrm{Det1}:=\{{\hat{t}}_{j}=t_{j,*}+\alpha\}=\{\lambda_{\max}(\frac{1}{\alpha}\sum_{t=t_{j,*}+1}^{t_{j,*}+\alpha}(\bm{I}-\hat{\bm{P}}_{\ast}\hat{\bm{P}}_{\ast}{}^{\prime})\hat{\bm{\ell}}_{t}\hat{\bm{\ell}}_{t}^{\prime}(\bm{I}-\hat{\bm{P}}_{\ast}\hat{\bm{P}}_{\ast}{}^{\prime}))>\omega_{evals}\}$, 2. 2.

$\mathrm{Proj{SVD}}:=\cap_{k=1}^{K}\mathrm{Proj{SVD}}_{k}$* where $\mathrm{Proj{SVD}}_{k}:=\{\mathrm{SE}([\hat{\bm{P}}_{\ast},\hat{\bm{P}}_{{\mathrm{rot}},k}])\leq\zeta_{k}^{+}\}$,* 3. 3.

$\mathrm{Del}:=\{\mathrm{SE}(\hat{\bm{P}},\bm{P})\leq\tilde{\varepsilon}\}$, 4. 4.

\mathrm{NoFalseDets}:=\{\text{for all \mathcal{J}^{\alpha}\subseteq[{\hat{t}}{j,fin},t{j+1}), }\\ \lambda_{\max}(\frac{1}{\alpha}\sum_{t\in\mathcal{J}^{\alpha}}(\bm{I}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime})\hat{\bm{\ell}}_{t}\hat{\bm{\ell}}_{t}^{\prime}(\bm{I}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime}))\leq\omega_{evals}\}** 5. 5.

$\Gamma_{0,\mathrm{end}}:=\{\mathrm{SE}(\hat{\bm{P}}_{*},\bm{P}_{*})\leq\tilde{\varepsilon}\}$, 6. 6.

\Gamma_{j,\mathrm{end}}:=\Gamma_{j-1,\mathrm{end}}\cap\big{(}(\mathrm{Det0}\cap\mathrm{Proj{SVD}}\cap\mathrm{Del}\cap\mathrm{NoFalseDets})\cup(\overline{\mathrm{Det0}}\cap\mathrm{Det1}\cap\mathrm{Proj{SVD}}\cap\mathrm{Del}\cap\mathrm{NoFalseDets})\big{)}.

Let $p_{0}$ denote the probability that, conditioned on $\Gamma_{j-1,\mathrm{end}}$, the change got detected at $t=t_{j,*}$, i.e., let

[TABLE]

Thus, $\Pr(\overline{\mathrm{Det0}}|\Gamma_{j-1,\mathrm{end}})=1-p_{0}$. It is not easy to bound $p_{0}$. However, as we will see, this will not be needed.

Assume that $\Gamma_{j-1,\mathrm{end}}\cap\overline{\mathrm{Det0}}$ holds. Consider the interval $\mathcal{J}^{\alpha}:=[t_{j,*},t_{j,*}+\alpha)$. This interval starts at or after $t_{j}$, so, for all $t$ in this interval, the subspace has changed. For this interval, $\bm{\Psi}=\bm{\Phi}=\bm{I}-\hat{\bm{P}}_{\ast}\hat{\bm{P}}_{\ast}{}^{\prime}$. Applying the last item of Theorem 5.7, w.p. at least $1-12n^{-12}$,

[TABLE]

where $q_{\mathrm{rot}}$ is the bound $\|(\bm{\Psi}_{\hat{\mathcal{T}}_{t}}{}^{\prime}\bm{\Psi}_{\hat{\mathcal{T}}_{t}})^{-1}\bm{I}_{\hat{\mathcal{T}}_{t}}{}^{\prime}\bm{\Psi}\bm{P}_{\mathrm{rot}}\|$. Theorem 5.7 is applicable for the reasons given in the proof of Lemma 5.16. Proceeding as in the proof of Lemma 5.16 for $k=1$, we get that $q_{\mathrm{rot}}=1.2((0.1+\tilde{\varepsilon})|\sin\theta|+\tilde{\varepsilon})$. Thus, using the bound on $\tilde{\varepsilon}$, we can conclude that, w.p. at least $1-12n^{-12}$,

[TABLE]

and thus ${\hat{t}}_{j}=t_{j,*}+\alpha$. This follows since $\omega_{evals}=5\tilde{\varepsilon}^{2}\lambda^{+}=5\tilde{\varepsilon}^{2}f\lambda^{-}\leq 5\tilde{\varepsilon}^{2}f^{2}\lambda^{-}\leq 5(0.01\min_{j}\mathrm{SE}(\bm{P}_{j-1},\bm{P}_{j}))^{2}\lambda^{-}$ and $\sin\theta=\mathrm{SE}(\bm{P}_{j-1},\bm{P}_{j})$. In other words,

[TABLE]

Conditioned on $\Gamma_{j-1,\mathrm{end}}\cap\overline{\mathrm{Det0}}\cap\mathrm{Det1}$, the first projection-SVD step is done at $t={\hat{t}}_{j}+\alpha=t_{j,*}+2\alpha$ and so on. We can state and prove Lemma 5.16 with $\Gamma_{k}$ replaced by $\Gamma_{j,\mathrm{end}}\cap\overline{\mathrm{Det0}}\cap\mathrm{Det1}\cap\mathrm{Proj{SVD}}_{1}\cap\mathrm{Proj{SVD}}_{2}\dots\mathrm{Proj{SVD}}_{k}$ and with the $k$-th projection-SVD interval being $\mathcal{J}_{k}:=[{\hat{t}}_{j}+(k-1)\alpha,{\hat{t}}_{j}+k\alpha)$. We can state and prove a similarly changed version of Lemma 5.17 for the simple SVD based deletion step. Applying Lemma 5.16 for each $k$, and then apply Lemma 5.17,

[TABLE]

We can also do a similar thing for the case when the change is detected at $t_{j,*}$, i.e. when $\mathrm{Det0}$ holds. In this case, we replace $\Gamma_{k}$ by $\Gamma_{j,\mathrm{end}}\cap\mathrm{Det0}\cap\mathrm{Proj{SVD}}_{1}\cap\mathrm{Proj{SVD}}_{2}\dots\mathrm{Proj{SVD}}_{k}$ and conclude that

[TABLE]

Finally consider the $\mathrm{NoFalseDets}$ event. First, assume that $\Gamma_{j-1,\mathrm{end}}\cap\mathrm{Det0}\cap\mathrm{Proj{SVD}}\cap\mathrm{Del}$ holds. Consider any interval $\mathcal{J}^{\alpha}\subseteq[{\hat{t}}_{j,fin},t_{j+1})$. In this interval, $\hat{\bm{P}}_{(t)}=\hat{\bm{P}}$, $\bm{\Psi}=\bm{\Phi}=\bm{I}-\hat{\bm{P}}\hat{\bm{P}}{}^{\prime}$ and $\mathrm{SE}(\hat{\bm{P}},\bm{P})\leq\tilde{\varepsilon}$. Also, using Lemma 5.15, $\bm{e}_{t}$ satisfies (10) for $t$ in this interval. Thus, defining $\bm{e}_{l,t}=\bm{\Phi}\bm{I}_{\mathcal{T}_{t}}(\bm{\Psi}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}_{\mathcal{T}_{t}})^{-1}\bm{I}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}\bm{\ell}_{t}$, $\bm{e}_{v,t}=\bm{\Phi}\bm{I}_{\mathcal{T}_{t}}(\bm{\Psi}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}_{\mathcal{T}_{t}})^{-1}\bm{I}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}\bm{v}_{t}$ and $\bm{z}_{t}=\bm{e}_{v,t}+\bm{\Phi}\bm{v}_{t}$

[TABLE]

We can bound the first term using Vershynin’s sub-Gaussian result (Theorem E.28) and the other terms using matrix Bernstein (Theorem E.27). The approach is similar to that of the proof of Lemma B.24. The derivation is more straightforward in this case, since for the above interval $\|\bm{\Psi}\bm{P}\|=\|\bm{\Phi}\bm{P}\|\leq\tilde{\varepsilon}$. The required bounds on $\alpha$ are also the same as those needed for Lemma B.24 to hold. We conclude that, w.p. at least $1-12n^{-12}$,

[TABLE]

This follows since $\omega_{evals}=5\tilde{\varepsilon}^{2}\lambda^{+}=5\tilde{\varepsilon}^{2}f\lambda^{-}$. Since $\mathrm{Det0}$ holds, ${\hat{t}}_{j}=t_{j,*}$. Thus, we have a total of $\lfloor\frac{t_{j+1}-t_{j,*}-K\alpha-\alpha}{\alpha}\rfloor$ intervals $\mathcal{J}^{\alpha}$ that are subsets of $[{\hat{t}}_{j,fin},t_{j+1})$. Moreover, $\lfloor\frac{t_{j+1}-t_{j,*}-K\alpha-\alpha}{\alpha}\rfloor\leq\lfloor\frac{t_{j+1}-t_{j}-K\alpha-\alpha}{\alpha}\rfloor\leq\lfloor\frac{t_{j+1}-t_{j}}{\alpha}\rfloor-(K+1)$ since $\alpha\leq\alpha$. Thus,

[TABLE]

On the other hand, if we condition on $\Gamma_{j-1,\mathrm{end}}\cap\overline{\mathrm{Det0}}\cap\mathrm{Det1}\cap\mathrm{Proj{SVD}}\cap\mathrm{Del}$, then ${\hat{t}}_{j}=t_{j,*}+\alpha$. Thus,

[TABLE]

We can now combine the above facts to bound $\Pr(\Gamma_{j,\mathrm{end}}|\Gamma_{j-1,\mathrm{end}})$. Recall that $p_{0}:=\Pr(\mathrm{Det0}|\Gamma_{j-1,\mathrm{end}})$. Clearly, the events $(\mathrm{Det0}\cap\mathrm{Proj{SVD}}\cap\mathrm{Del}\cap\mathrm{NoFalseDets})$ and $(\overline{\mathrm{Det0}}\cap\mathrm{Det1}\cap\mathrm{Proj{SVD}}\cap\mathrm{Del}\cap\mathrm{NoFalseDets})$ are disjoint. Thus,

[TABLE]

Since the events $\Gamma_{j,\mathrm{end}}$ are nested, the above implies that

[TABLE]

Appendix B Proof of Theorem 5.7: PCA in data-dependent noise with partial subspace knowledge

We prove Theorem 5.7 with the modification given in Remark 5.8. Thus we condition on $\mathcal{E}_{0}$ defined in the remark. Recall that $\bm{\Phi}:=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}$. Let

[TABLE]

denote the reduced QR decomposition of $(\bm{\Phi}\bm{P}_{{\mathrm{rot}}})$. Here, and in the rest of this proof, we write things in a general fashion to allow $\bm{P}_{\mathrm{rot}}$ to contain more than one direction. This makes it easier to understand how our guarantees extend to the more general case ($\bm{P}_{\mathrm{rot}}$ being an $n\times r_{\mathrm{ch}}$ basis matrix with $r_{\mathrm{ch}}>1$) easier. The proof uses the following simple lemma at various places.

Lemma B.21**.**

Assume that $\mathcal{E}_{0}$ holds. Then,

$\left\|\bm{M}_{1,t}\bm{P}_{{\mathrm{fix}}}\right\|\leq q_{0}$*, $\left\|\bm{M}_{1,t}\bm{P}_{{\mathrm{ch}}}\right\|\leq q_{0}$ and $\left\|\bm{M}_{1,t}\bm{P}_{{\mathrm{rot}}}\right\|\leq q_{{\mathrm{rot}}}$ * 2. 2.

$\left\|\bm{\Phi}\bm{P}_{{\mathrm{fix}}}\right\|\leq\tilde{\varepsilon}$, $\left\|\bm{\Phi}\bm{P}_{{\mathrm{ch}}}\right\|\leq\tilde{\varepsilon}$, $\left\|\bm{\Phi}\bm{P}_{\mathrm{new}}\right\|\leq 1$, 3. 3.

$\left\|\bm{R}_{{\mathrm{rot}}}\right\|=\left\|\bm{\Phi}\bm{P}_{{\mathrm{rot}}}\right\|\leq\tilde{\varepsilon}|\cos\theta|+|\sin\theta|\leq\tilde{\varepsilon}+|\sin\theta|$** 4. 4.

$\sigma_{\min}(\bm{R}_{{\mathrm{rot}}})=\sigma_{\min}(\bm{\Phi}\bm{P}_{{\mathrm{rot}}})\geq\sqrt{\sin^{2}\theta(1-\tilde{\varepsilon}^{2})-2\tilde{\varepsilon}|\sin\theta|}$** 5. 5.

$\left\|\bm{\Phi}\bm{\ell}_{t}\right\|\leq 2\tilde{\varepsilon}\sqrt{\eta r\lambda^{+}}+|\sin\theta|\sqrt{\eta\lambda_{{\mathrm{ch}}}}$.

Proof of Lemma B.21.

Item 1 follows because $\left\|\bm{M}_{1,t}\bm{P}_{*}\right\|\leq q_{0}\leq 2\tilde{\varepsilon}$ and $\left\|\bm{M}_{1,t}\bm{P}_{*}\right\|=\left\|\bm{M}_{1,t}[\bm{P}_{{\mathrm{fix}}},\bm{P}_{{\mathrm{ch}}}]\right\|\geq\left\|\bm{M}_{1,t}\bm{P}_{{\mathrm{fix}}}\right\|$. Similarly, $\left\|\bm{M}_{1,t}\bm{P}_{*}\right\|\geq\left\|\bm{M}_{1,t}\bm{P}_{{\mathrm{ch}}}\right\|$. The first two claims of item 2 follow because $\left\|\bm{\Phi}\bm{P}_{*}\right\|\leq\tilde{\varepsilon}$ and the bound on item 1. Third claim uses $\|\bm{\Phi}\bm{P}_{\mathrm{new}}\|\leq\|\bm{\Phi}\|\|\bm{P}_{\mathrm{new}}\|=1$. The fourth claim uses triangle inequality and definition of $\bm{P}_{\mathrm{rot}}$. For Item 3, recall that $\bm{\Phi}\bm{P}_{{\mathrm{rot}}}\overset{\mathrm{QR}}{=}\bm{E}_{{\mathrm{rot}}}\bm{R}_{{\mathrm{rot}}}$. Thus, $\sigma_{i}(\bm{R}_{{\mathrm{rot}}})=\sigma_{i}(\bm{\Phi}\bm{P}_{{\mathrm{rot}}})$. Thus $\|\bm{R}_{{\mathrm{rot}}}\|=\left\|\bm{\Phi}\bm{P}_{{\mathrm{rot}}}\right\|\leq\tilde{\varepsilon}+|\sin\theta|$

Item 4: From above, $\sigma_{\min}(\bm{R}_{{\mathrm{rot}}})=\sigma_{\min}(\bm{\Phi}\bm{P}_{{\mathrm{rot}}})$. Moreover, $\sigma_{\min}(\bm{\Phi}\bm{P}_{{\mathrm{rot}}})=\sqrt{\lambda_{\min}(\bm{P}_{{\mathrm{rot}}}{}^{\prime}\bm{\Phi}^{\prime}\bm{\Phi}\bm{P}_{{\mathrm{rot}}})}=\sqrt{\lambda_{\min}(\bm{P}_{{\mathrm{rot}}}{}^{\prime}\bm{\Phi}\bm{P}_{{\mathrm{rot}}})}$. We bound this as follows. Recall that $\bm{\Phi}=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}$.

[TABLE]

The last inequality used Lemma 5.10.

Item 5: Using the previous items and the definition of $\eta$,

[TABLE]

∎

B-A Proof of Theorem 5.7

Proof of Theorem 5.7.

We have

[TABLE]

where the last inequality used Lemma B.21. Consider $\mathrm{SE}([\hat{\bm{P}}_{*},\hat{\bm{P}}_{\mathrm{rot}}],\bm{P}_{{\mathrm{rot}}})$.

[TABLE]

The last inequality used Lemma B.21. To bound $\mathrm{SE}(\hat{\bm{P}}_{\mathrm{rot}},\bm{E}_{{\mathrm{rot}}})$, we use the Davis-Kahan $\sin\theta$ theorem [33] given below.

Theorem B.22** (Davis-Kahan $\sin\theta$ theorem).**

Consider $n\times n$ Hermitian matrices, $\bm{D}$ and $\hat{\bm{D}}$ such that

[TABLE]

where $[\bm{E},\bm{E}_{\perp}]$ and $[\bm{F},\bm{F}_{\perp}]$ are orthogonal matrices and $\mathrm{rank}(\bm{F})=\mathrm{rank}(\bm{E})$. Let

[TABLE]

If $\lambda_{\min}(\bm{A})-\lambda_{\max}(\bm{A}_{\mathrm{rest}})-\|\bm{H}\|>0$ and $\mathrm{rank}(\bm{E})=\mathrm{rank}(\bm{F})$, then,

[TABLE]

To use this result to bound $\mathrm{SE}(\hat{\bm{P}}_{\mathrm{rot}},\bm{E}_{\mathrm{rot}})$, let $\hat{\bm{D}}:=\bm{D}_{obs}=\frac{1}{\alpha}\sum_{t}\bm{\Phi}\bm{y}_{t}\bm{y}_{t}{}^{\prime}\bm{\Phi}$. Its top eigenvector is $\hat{\bm{P}}_{\mathrm{rot}}$. We need to define a matrix $\bm{D}$ that is such that its top eigenvector is $\bm{E}_{{\mathrm{rot}}}$ and the gap between its first and second eigenvalues is more than $\|\bm{H}\|$. Consider the matrix

[TABLE]

If $\lambda_{\max}(\bm{A}_{\mathrm{rest}})<\lambda_{\min}(\bm{A})$, then $\bm{E}_{{\mathrm{rot}}}$ is the top eigenvector of $\bm{D}$. Moreover, if $\lambda_{\max}(\bm{A}_{\mathrm{rest}})<\lambda_{\min}(\bm{A})-\|\bm{H}\|$, then the gap requirement holds too. Thus, by the $\sin\theta$ theorem,

[TABLE]

Here again, we should point out that, in the simple case that we consider where $\bm{P}_{\mathrm{rot}}$ is a vector (only one direction changes), $\bm{A}$ is a non-negative scalar and $\lambda_{\min}(\bm{A})=\bm{A}$. However the above discussion applies even in the general case when $r_{\mathrm{ch}}>1$. The rest of the proof obtains high probability bounds on the terms in the above expression. $\left\|\bm{H}\right\|=\left\|\bm{D}-\hat{\bm{D}}\right\|$ can be bounded as follows.

Lemma B.23**.**

*Let

${\mathrm{term11}}=\frac{1}{\alpha}\sum_{t}\bm{E}_{{\mathrm{rot}}}\bm{E}_{{\mathrm{rot}}}{}^{\prime}\bm{\Phi}\bm{\ell}_{t}\bm{\ell}_{t}{}^{\prime}\bm{\Phi}\bm{E}_{{\mathrm{rot}},\perp}\bm{E}_{{\mathrm{rot}},\perp}{}^{\prime}$. Then,*

[TABLE]

Proof of Lemma B.23.

Recall that $\bm{H}=\hat{\bm{D}}-\bm{D}$. Thus

[TABLE]

Using triangle inequality the bound follows. ∎

The next lemma obtains high probability bounds on the above terms and the two other terms from (19).

Lemma B.24**.**

Assume that the assumptions of Theorem 5.7 with the modification given in Remark 5.8 hold. Let $\epsilon_{0}=0.01|\sin\theta|(\tilde{\varepsilon}+q_{{\mathrm{rot}}})$, $\epsilon_{1}=0.01(q_{{\mathrm{rot}}}^{2}+\tilde{\varepsilon}^{2})$, and $\epsilon_{2}=0.01$. For an $\alpha\geq\alpha_{0}:=C\eta\max\left\{fr\log n,\ \eta f^{2}(r+\log n)\right\}$, conditioned on $\mathcal{E}_{0}$, all the following hold w.p. at least $1-12n^{-12}$:

$\left\|\frac{1}{\alpha}\sum_{t}\bm{\Phi}\bm{\ell}_{t}\bm{w}_{t}{}^{\prime}\bm{\Phi}\right\|\leq\left[\sqrt{b_{0}}\left(2\tilde{\varepsilon}^{2}f+(\tilde{\varepsilon}+|\sin\theta|)q_{{\mathrm{rot}}}\right)+\epsilon_{0}\right]\lambda_{{\mathrm{ch}}}$, 2. 2.

$\left\|\frac{1}{\alpha}\sum_{t}\bm{\Phi}\bm{w}_{t}\bm{w}_{t}{}^{\prime}\bm{\Phi}\right\|\leq\left[\sqrt{b_{0}}\left(4\tilde{\varepsilon}^{2}f+q_{{\mathrm{rot}}}^{2}\right)+\epsilon_{1}\right]\lambda_{{\mathrm{ch}}}$, 3. 3.

$\lambda_{\min}(\bm{A})\geq(\sin^{2}\theta(1-\tilde{\varepsilon}^{2})-2\tilde{\varepsilon}|\sin\theta|)(1-\epsilon_{2})\lambda_{{\mathrm{ch}}}-2\tilde{\varepsilon}(\tilde{\varepsilon}+|\sin\theta|)\epsilon_{2}\lambda_{{\mathrm{ch}}}$, 4. 4.

$\lambda_{\max}(\bm{A}_{\mathrm{rest}})\leq\tilde{\varepsilon}^{2}\lambda^{+}+\tilde{\varepsilon}^{2}\epsilon_{2}\lambda_{{\mathrm{ch}}}$, 5. 5.

$\left\|{\mathrm{term11}}\right\|\leq\left[\tilde{\varepsilon}^{2}f+2\tilde{\varepsilon}^{2}\epsilon_{2}+\tilde{\varepsilon}|\sin\theta|\epsilon_{2}|\right]\lambda_{{\mathrm{ch}}}.$** 6. 6.

$\left\|\frac{1}{\alpha}\sum_{t}\bm{\Phi}\bm{\ell}_{t}\bm{z}_{t}{}^{\prime}\bm{\Phi}\right\|\leq\epsilon_{0}\lambda_{{\mathrm{ch}}}$** 7. 7.

$\left\|\frac{1}{\alpha}\sum_{t}\bm{\Phi}\bm{z}_{t}\bm{z}_{t}{}^{\prime}\bm{\Phi}\right\|\leq\left[\left(4\tilde{\varepsilon}^{2}f+q_{{\mathrm{rot}}}^{2}\right)+\epsilon_{1}\right]\lambda_{{\mathrm{ch}}}$**

Using Lemma B.24 and substituting for $\epsilon_{0},\epsilon_{1},\epsilon_{2}$, we conclude the following. Conditioned on $\mathcal{E}_{0}$, with probability at least $1-12n^{-12}$,

[TABLE]

where $\mathrm{numer}$ denotes the numerator expression. The numerator, $\mathrm{numer}$, expression can be simplified to

[TABLE]

Further, using $\tilde{\varepsilon}f\leq 0.01|\sin\theta|$, $\sqrt{b_{0}}\leq 0.1$ and $q_{{\mathrm{rot}}}\leq 0.2|\sin\theta|$,

[TABLE]

This can be loosely upper bounded by $0.26\sin^{2}\theta$. We use this loose upper bound when this term appears in the denominator. Following a similar approach for the denominator, denoted $\mathrm{denom}$,

[TABLE]

Thus,

[TABLE]

Using (17) and $\tilde{\varepsilon}\leq\tilde{\varepsilon}f\leq 0.01|\sin\theta|$,

[TABLE]

*Proof of the last claim: lower bound on $\lambda_{\max}(\bm{D}_{obs})$. * Using Weyl’s inequality,

[TABLE]

Using the bounds from Lemmas B.23 and B.24 and (12), we get the lower bound. ∎

B-B Proof of Lemma B.24: high probability bounds on the $\sin\theta$ theorem bound terms

Proof of Lemma B.24.

Recall the definition of the event $\mathcal{E}_{0}$ from Remark 5.8. To prove this lemma, we first bound the probabilities of all the events conditioned on $\{\hat{\bm{P}}_{*},Z\}$, for values of $\{\hat{\bm{P}}_{*},Z\}\in\mathcal{E}_{0}$. Then we use the following simple fact.

Fact B.25**.**

If $\Pr(\mathrm{Event}|\{\hat{\bm{P}}_{*},Z\})\geq p_{0}$ for all $\{\hat{\bm{P}}_{*},Z\}\in\mathcal{E}_{0}$, then,

[TABLE]

In the discussion below, we condition on $\{\hat{\bm{P}}_{*},Z\}$, for values of $\{\hat{\bm{P}}_{*},Z\}$ in $\mathcal{E}_{0}$. Conditioned on $\{\hat{\bm{P}}_{*},Z\}$, the matrices $\bm{E}_{\mathrm{rot}}$, $\bm{E}_{{\mathrm{rot}},\perp}$, $\bm{\Phi}$, etc, are constants (not random). All the terms that we bound in this lemma are either of the form $\sum_{t\in\mathcal{J}^{\alpha}}g_{1}(\hat{\bm{P}}_{*},Z)\bm{\ell}_{t}\bm{\ell}_{t}{}^{\prime}g_{2}(\hat{\bm{P}}_{*},Z)$, for some functions $g_{1}(.),g_{2}(.)$, or are sub-matrices of such a term.

Since the pair $\{\hat{\bm{P}}_{*},Z\}$ is independent of the $\bm{\ell}_{t}$’s for $t\in\mathcal{J}^{\alpha}$, and these $\bm{\ell}_{t}$’s are mutually independent, hence, even conditioned on $\{\hat{\bm{P}}_{*},Z\}$, the same holds: the $\bm{\ell}_{t}$’s for $t\in\mathcal{J}^{\alpha}$ are mutually independent. Thus, once we condition on $\{\hat{\bm{P}}_{*},Z\}$, the summands in the terms we need to bound are mutually independent. As a result, matrix Bernstein (Theorem E.27) or Vershynin’s sub-Gaussian result (Theorem E.28) are applicable.

Item 1: In the proof of this and later items, we condition on $\{\hat{\bm{P}}_{*},Z\}$, for values of $\{\hat{\bm{P}}_{*},Z\}$ in $\mathcal{E}_{0}$.

Since $\left\|\bm{\Phi}\right\|=1$,

[TABLE]

To bound the RHS above, we will apply matrix Bernstein (Theorem E.27) with $\bm{Z}_{t}=\bm{\Phi}\bm{\ell}_{t}\bm{w}_{t}{}^{\prime}$. As explained above, conditioned on $\{\hat{\bm{P}}_{*},Z\}$, the $\bm{Z}_{t}$’s are mutually independent. We first obtain a bound on the expected value of the time average of the $\bm{Z}_{t}$’s and then compute $R$ and $\sigma^{2}$ needed by Theorem E.27. By Cauchy-Schwarz,

[TABLE]

where (a) follows by Cauchy-Schwarz (Theorem E.26) with $\bm{X}_{t}=\bm{\Phi}\bm{P}\bm{\Lambda}\bm{P}{}^{\prime}\bm{M}_{1,t}{}^{\prime}$ and $\bm{Y}_{t}=\bm{M}_{2,t}$, (b) follows from the assumption on $\bm{M}_{2,t}$, and the last inequality follows from Lemma B.21. Using $q_{0}\leq 2\tilde{\varepsilon}$,

[TABLE]

To compute $R$, using Lemma B.21 and using $q_{0}\leq 2\tilde{\varepsilon}$ and $q_{\mathrm{rot}}<|\sin\theta|$,

[TABLE]

for numerical constants $c_{1},c_{2}$. Next we compute $\sigma^{2}$. Since $\bm{w}_{t}$’s are bounded r.v.’s, we have

[TABLE]

for numerical constants $c_{1}$ and $c_{2}$. The above bounds again used $q_{0}\leq 2\tilde{\varepsilon}$ and $q_{\mathrm{rot}}<|\sin\theta|$. For bounding $\left\|\frac{1}{\alpha}\sum_{t}\mathbb{E}{\left[\bm{Z}_{t}{}^{\prime}\bm{Z}_{t}\right]}\right\|$ we get the same expression except for the values of $c_{1}$, $c_{2}$. Thus, applying matrix Bernstein (Theorem E.27) followed by Fact B.25,

[TABLE]

Let $\epsilon=\epsilon_{0}\lambda_{{\mathrm{ch}}}$ where $\epsilon_{0}=0.01\sin\theta(q_{{\mathrm{rot}}}+\tilde{\varepsilon})$. Then, clearly,

[TABLE]

Hence, for the probability to be of the form $1-2n^{-12}$ we require that $\alpha\geq\alpha_{(1)}$ where

[TABLE]

Thus, if $\alpha\geq\alpha_{(1)}$, conditioned on $\mathcal{E}_{0}$, the bound on $\left\|\frac{1}{\alpha}\sum_{t}\bm{\Phi}\bm{\ell}_{t}\bm{w}_{t}{}^{\prime}\bm{\Phi}\right\|$ given in Lemma B.24 holds w.p. at least $1-2n^{-12}$.

Item 2: We use Theorem E.27 (matrix Bernstein) with $\bm{Z}_{t}:=\bm{\Phi}\bm{w}_{t}\bm{w}_{t}{}^{\prime}\bm{\Phi}$. The proof approach is similar to that of the proof of item 1. First we bound the norm of the expectation of the time average of $\bm{Z}_{t}$:

[TABLE]

(a) follows from Cauchy-Schwarz (Theorem E.26) with $\bm{X}_{t}=\bm{M}_{2,t}$ and $\bm{Y}_{t}=\bm{M}_{1,t}\bm{P}\bm{\Lambda}\bm{P}{}^{\prime}\bm{M}_{1,t}{}^{\prime}\bm{M}_{2,t}{}^{\prime}$, (b) follows from the assumption on $\bm{M}_{2,t}$, and (c) follows from Lemma B.21. The last inequality used $q_{0}\leq 2\tilde{\varepsilon}$. To obtain $R$,

[TABLE]

To obtain $\sigma^{2}$,

[TABLE]

Applying matrix Bernstein (Theorem E.27) followed by Fact B.25, we have

[TABLE]

Let $\epsilon=\epsilon_{1}\lambda_{{\mathrm{ch}}}$, $\epsilon_{1}=0.01(q_{{\mathrm{rot}}}^{2}+\tilde{\varepsilon}^{2})$. Then we get

[TABLE]

For the success probability to be of the form $1-2n^{-12}$ we require $\alpha\geq\alpha_{(2)}$ where

[TABLE]

Thus, if $\alpha\geq\alpha_{(2)}$,

$\Pr\left(\left\|\frac{1}{\alpha}\sum_{t}\bm{\Phi}\bm{w}_{t}\bm{w}_{t}{}^{\prime}\bm{\Phi}\right\|\leq\left[\sqrt{b_{0}}\left(4\tilde{\varepsilon}^{2}f+q_{{\mathrm{rot}}}^{2}\right)+\epsilon_{1}\right]\lambda_{{\mathrm{ch}}}|\mathcal{E}_{0}\right)\geq 1-n^{-12}.$

Item 3: Expanding the expression for $\bm{A}$,

[TABLE]

where $\mathrm{term1}:=\bm{E}_{{\mathrm{rot}}}{}^{\prime}\bm{\Phi}\bm{P}_{{\mathrm{fix}}}\left(\frac{1}{\alpha}\sum_{t}\bm{a}_{t,{\mathrm{fix}}}\bm{a}_{t,{\mathrm{ch}}}{}^{\prime}\right)\bm{P}_{{\mathrm{rot}}}{}^{\prime}\bm{\Phi}\bm{E}_{{\mathrm{rot}}}$. Since the first term on the RHS is positive semi-definite,

[TABLE]

Under our current assumptions, the $\bm{a}_{t,{\mathrm{ch}}}$’s are scalars, so $\bm{A}$ and $\frac{1}{\alpha}\sum_{t}\bm{a}_{t,{\mathrm{ch}}}\bm{a}_{t,{\mathrm{ch}}}{}^{\prime}$ are actually scalars. However, we write things in a general fashion (allowing $\bm{a}_{t,{\mathrm{ch}}}$’s to be $r_{\mathrm{ch}}$ length vectors), so as to make our later discussion of the $r_{\mathrm{ch}}>1$ case easier. Using (21),

[TABLE]

The second inequality follows using $\bm{E}_{{\mathrm{rot}}}{}^{\prime}\bm{\Phi}\bm{P}_{{\mathrm{rot}}}=\bm{E}_{{\mathrm{rot}}}{}^{\prime}\bm{E}_{{\mathrm{rot}}}\bm{R}_{\mathrm{rot}}=\bm{R}_{\mathrm{rot}}$ and Lemma B.21. The first inequality is straightforward if $\bm{a}_{t,{\mathrm{ch}}}$’s are scalars (current setting); it follows using Ostrowski’s theorem [35] in the general case.

To bound the remaining terms in the above expression, we use Vershynin’s sub-Gaussian result [32, Theorem 5.39] summarized in Theorem E.28. To apply this, recall that $(\bm{a}_{t})_{i}$ are bounded random variables satisfying $|(\bm{a}_{t})_{i}|\leq\sqrt{\eta\lambda_{i}}$. Hence they are sub-Gaussian with sub-Gaussian norm $\sqrt{\eta\lambda_{i}}$ [32]. Using [32, Lemma 5.24], the vectors $\bm{a}_{t}$ are also sub-Gaussian with sub-Gaussian norm bounded by $\max_{i}\sqrt{\eta\lambda_{i}}=\sqrt{\eta\lambda^{+}}$. Thus, applying Theorem E.28 with $K\equiv\sqrt{\eta\lambda^{+}}$, $\epsilon\equiv\epsilon_{2}\lambda_{\mathrm{ch}}$, $N\equiv\alpha$, $n_{w}\equiv r$, followed by using Fact B.25, if $\alpha\geq\alpha_{(3)}:=\frac{C(r\log 9+10\log n)f^{2}}{\epsilon_{2}^{2}},$ then,

[TABLE]

We could also have used matrix Bernstein to bound $\left\|\sum_{t}\bm{a}_{t}\bm{a}_{t}{}^{\prime}\right\|$. However, since the $\bm{a}_{t}$’s are $r$-length vectors and $r\ll n$, the Vershynin result requires a smaller lower bound on $\alpha$.

If $\bm{B}_{1}$ is a sub-matrix of a matrix $\bm{B}$, then $\left\|\bm{B}_{1}\right\|\leq\left\|\bm{B}\right\|$. Thus, we can also use (23) for bounding the norm of various sub-matrices of $\left(\frac{1}{\alpha}\sum_{t}\bm{a}_{t}\bm{a}_{t}{}^{\prime}-\bm{\Lambda}\right)$. Doing this, we get

[TABLE]

Combining (B-B), (B-B) and (B-B), if $\alpha\geq\alpha_{(3)}$,

[TABLE]

Item 4: Recall that $\bm{E}_{{\mathrm{rot}},\perp}{}^{\prime}\bm{\Phi}\bm{P}_{{\mathrm{rot}}}=0$. Thus,

[TABLE]

where the last inequality follows from Ostrowski’s theorem [35]. Using this and (B-B), if $\alpha\geq\alpha_{(3)}$,

[TABLE]

Item 5: Recall that

${\mathrm{term11}}~{}=~{}\frac{1}{\alpha}\sum\bm{E}_{{\mathrm{rot}}}\bm{E}_{{\mathrm{rot}}}{}^{\prime}\bm{\Phi}\bm{\ell}_{t}\bm{\ell}_{t}{}^{\prime}\bm{\Phi}\bm{E}_{{\mathrm{rot}},\perp}\bm{E}_{{\mathrm{rot}},\perp}{}^{\prime}$. As in earlier items, we can expand this into a sum of four terms using $\bm{\ell}_{t}=\bm{P}_{{\mathrm{fix}}}\bm{a}_{t,{\mathrm{fix}}}+\bm{P}_{{\mathrm{rot}}}\bm{a}_{t,{\mathrm{ch}}}$. Then using $\bm{E}_{{\mathrm{rot}},\perp}{}^{\prime}\bm{\Phi}\bm{P}_{{\mathrm{rot}}}=0$ and $\left\|\bm{E}_{{\mathrm{rot}}}\right\|=\left\|\bm{E}_{{\mathrm{rot}},\perp}\right\|=1$, we get

[TABLE]

Using (B-B) and (B-B), if $\alpha\geq\alpha_{(3)}$, w.p. at least $1-4n^{-12}$, conditioned on $\mathcal{E}_{0}$, $\left\|{\mathrm{term11}}\right\|\leq\tilde{\varepsilon}^{2}(\lambda^{+}+\epsilon_{2}\lambda_{{\mathrm{ch}}})+(\tilde{\varepsilon}(\tilde{\varepsilon}+|\sin\theta|))\epsilon_{2}\lambda_{{\mathrm{ch}}}.$

Item 6: Consider $\left\|\frac{1}{\alpha}\sum_{t}\bm{\Phi}\bm{\ell}_{t}\bm{z}_{t}{}^{\prime}\right\|$. We will apply matrix Bernstein (Theorem E.27). We have $\left\|\mathbb{E}{\left[\frac{1}{\alpha}\sum_{t}\bm{\Phi}\bm{\ell}_{t}\bm{z}_{t}{}^{\prime}\right]}\right\|=0$ since $\bm{\ell}_{t}$’s are independent of $\bm{z}_{t}$’s and both are zero mean. We obtain $R$ as follows

[TABLE]

for numerical constants $c_{1},c_{2}$. Next we compute $\sigma^{2}$ as follows. First consider

[TABLE]

we note here that since $b_{z}^{2}=r\lambda_{z}^{+}$, the other term in the expression for $\sigma^{2}$ is the same (modulo constants) as $\sigma_{1}^{2}$. Furthermore, notice that the expressions for both $R$ and $\sigma^{2}$ are the same as the ones obtained in Item 1. Thus, we use the same deviation, $\epsilon_{0}$ here, and hence also obtain the same sample complexity, $\alpha$; i.e., we let $\epsilon=\epsilon_{0}\lambda_{{\mathrm{ch}}}$ where $\epsilon_{0}=0.01\sin\theta(q_{{\mathrm{rot}}}+\tilde{\varepsilon})$, and obtain $\alpha\geq\alpha_{(1)}$ derived in item 1.

Item 7: This term follows in a similar fashion as Item 2, 6 and $\alpha\geq\alpha_{(2)}$ suffices. ∎

Appendix C Proof of Projected CS Lemma

Proof of Lemma 5.15.

The first four claims were already proved below the lemma statement. Consider the fifth claim (exact support recovery). Recall that for any $t\in\mathcal{J}_{k}$, $\bm{v}_{t}$ satisfies $\left\|\bm{v}_{t}\right\|\leq C(2\tilde{\varepsilon}\sqrt{r\lambda^{+}}+\zeta_{k-1}^{+}\sqrt{\lambda_{{\mathrm{ch}}}}):=b_{v,t}$ (for $t\in\mathcal{J}_{1}$ and $t\in\mathcal{J}_{0}$ the bounds are the same) with $C=\sqrt{\eta}$, and thus $\bm{b}_{t}:=\bm{\Psi}(\bm{\ell}_{t}+\bm{v}_{t})$ satisfies

[TABLE]

From the lower bound on $x_{\min,t}$ in Theorem 2.2 or that in Corollary 2.3, $b_{b,t}<x_{\min,t}/15$. Also, we set $\xi_{t}=x_{\min,t}/15$. Using these facts, and $\delta_{2s}(\bm{\Psi})\leq 0.12<0.15$ (third claim of this lemma), [10, Theorem 1.2] implies that

[TABLE]

Thus,

[TABLE]

We have $\omega_{supp,t}=x_{\min,t}/2$. Consider an index $i\in\mathcal{T}_{t}$. Since $|(\bm{x}_{t})_{i}|\geq x_{\min,t}$,

[TABLE]

Thus, $|(\bm{\hat{x}}_{t,cs})_{i}|>\frac{x_{\min,t}}{2}=\omega_{supp,t}$ which means $i\in\hat{\mathcal{T}}_{t}$. Hence $\mathcal{T}_{t}\subseteq\hat{\mathcal{T}}_{t}$. Next, consider any $j\notin\mathcal{T}_{t}$. Then, $(\bm{x}_{t})_{j}=0$ and so

[TABLE]

which implies $j\notin\hat{\mathcal{T}}_{t}$ and so $\hat{\mathcal{T}}_{t}\subseteq\mathcal{T}_{t}$. Thus $\hat{\mathcal{T}}_{t}=\mathcal{T}_{t}$.

With $\hat{\mathcal{T}}_{t}=\mathcal{T}_{t}$, the sixth claim follows easily. Since $\mathcal{T}_{t}$ is the support of $\bm{x}_{t}$, $\bm{x}_{t}=\bm{I}_{\mathcal{T}_{t}}\bm{I}_{\mathcal{T}_{t}}{}^{\prime}\bm{x}_{t}$, and so

[TABLE]

since $\bm{\Psi}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}=\bm{I}_{\mathcal{T}_{t}}^{\prime}\bm{\Psi}^{\prime}\bm{\Psi}=\bm{I}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}$. Thus $\bm{e}_{t}=\hat{\bm{x}}_{t}-\bm{x}_{t}$ satisfies (10). Using (10) and the earlier claims,

[TABLE]

When $k=1$, $\bm{\Psi}=\bm{I}-\hat{\bm{P}}_{*}\hat{\bm{P}}_{*}{}^{\prime}$. Thus, using (13) and $\|\hat{\bm{P}}_{*}{}^{\prime}\bm{P}_{\mathrm{new}}\|\leq\tilde{\varepsilon}$ (follows from Lemma 5.10),

[TABLE]

also, in this interval, $b_{v,t}\leq 2\tilde{\varepsilon}\sqrt{\eta r\lambda^{+}}+0.11|\sin\theta|\sqrt{\eta\lambda_{\mathrm{ch}}}$ so that $\left\|\bm{e}_{t}\right\|\leq 2.4\cdot\left(2\tilde{\varepsilon}\sqrt{\eta r\lambda^{+}}+0.11|\sin\theta|\sqrt{\eta\lambda_{\mathrm{ch}}}\right)$ When $k>1$, $\left\|\bm{I}_{\mathcal{T}_{t}}{}^{\prime}\bm{\Psi}\bm{\ell}_{t}\right\|\leq\left\|\bm{\Psi}\bm{\ell}_{t}\right\|\leq\tilde{\varepsilon}\sqrt{r\eta\lambda^{+}}+\zeta_{k-1}^{+}\sqrt{\eta\lambda_{{\mathrm{ch}}}}.$ and the same bound holds on $b_{v,t}$ so that $\left\|\bm{e}_{t}\right\|\leq 2.4\cdot\left(\tilde{\varepsilon}\sqrt{r\eta\lambda^{+}}+\zeta_{k-1}^{+}\sqrt{\eta\lambda_{{\mathrm{ch}}}}\right)$ ∎

Appendix D Time complexity of s-ReProCS

The time-consuming steps of s-ReProCS are either $l_{1}$ minimization or the subspace update steps. Support estimation and LS steps are much faster and hence can be ignored for this discussion. The computational complexity of $l_{1}$ minimization (if the best solver were used) [36] is the cost of multiplying the CS matrix or its transpose with a vector times $\log(1/\epsilon)$ if $\epsilon$ is the bound on the error w.r.t. the true minimizer of the program. In ReProCS, the CS matrix is of the form $\bm{I}-\hat{\bm{P}}\hat{\bm{P}}^{\prime}$ where $\hat{\bm{P}}$ is of size $n\times r$ or $n\times(r+1)$, thus multiplying a vector with it takes time $O(nr)$. Thus, the $l_{1}$ minimization complexity per frame is $O(nr\log(1/\epsilon))$, and thus the total cost for $d-t_{\mathrm{train}}$ frames is $O(nr\log(1/\epsilon)(d-t_{\mathrm{train}}))$. The subspace update step consists of $(d-t_{\mathrm{train}}-J\alpha)/\alpha$ rank one SVD’s on an $n\times\alpha$ matrix (for either detecting subspace change or for projection-SVD) and $J$ rank $r$ SVD’s on an $n\times\alpha$ matrix (for subspace re-estimation). Thus the subspace update complexity is at most $O(n(d-t_{\mathrm{train}})r\log(1/\epsilon))$ and the total ReProCS complexity (without the initialization step) is $O(n(d-t_{\mathrm{train}})r\log(1/\epsilon))$.

If we assume that the initialization uses AltProj, AltProj is applied to a matrix of size $n\times t_{\mathrm{train}}$ with rank $r$. Thus the initialization complexity is $O(nt_{\mathrm{train}}r^{2}\log(1/\epsilon))$. If instead GD [8] is used, then the time complexity is reduced to $O(nt_{\mathrm{train}}rf\log(1/\epsilon))$. Treating $f$ as a constant (our discussion treats condition numbers as constants), the final complexity of s-ReProCS is $O(ndr\log(1/\epsilon))$.

If s-ReProCS is used to only solve the RPCA problem (compute column span of the entire matrix $\bm{L}$), then the SVD based subspace re-estimation step can be removed. With this change, the complexity of s-ReProCS (without the initialization step) reduces to just $O(nd\log(1/\epsilon))$ since only 1-SVDs are needed. Of course this would mean a slightly tighter bound on max-outlier-frac-col is required – it will need to be less than $c/(r+J)$.

Appendix E Preliminaries: Cauchy-Schwarz, matrix Bernstein and Vershynin’s sub-Gaussian result

Cauchy-Schwarz for sums of matrices says the following [5].

Theorem E.26**.**

For matrices $\bm{X}$ and $\bm{Y}$ we have

[TABLE]

Matrix Bernstein [31], conditioned on another r.v. $X$, says the following.

Theorem E.27**.**

Given an $\alpha$-length sequence of $n_{1}\times n_{2}$ dimensional random matrices and a r.v. $X$ Assume the following. For all $X\in\mathcal{C}$, (i) conditioned on $X$, the matrices $\bm{Z}_{t}$ are mutually independent, (i) $\mathbb{P}(\left\|\bm{Z}_{t}\right\|\leq R|X)=1$, and (iii) $\max\left\{\left\|\frac{1}{\alpha}\sum_{t}\mathbb{E}{\left[\bm{Z}_{t}{}^{\prime}\bm{Z}_{t}|X\right]}\right\|,\ \left\|\frac{1}{\alpha}\sum_{t}\mathbb{E}{\left[\bm{Z}_{t}\bm{Z}_{t}{}^{\prime}|X\right]}\right\|\right\}\leq\sigma^{2}$. Then, for an $\epsilon>0$,

[TABLE]

Vershynin’s result for matrices with independent sub-Gaussian rows [32, Theorem 5.39], conditioned on another r.v. $X$, says the following.

Theorem E.28**.**

Given an $N$-length sequence of sub-Gaussian random vectors $\bm{w}_{i}$ in $\mathbb{R}^{n_{w}}$, an r.v $X$, and a set $\mathcal{C}$. Assume that for all $X\in\mathcal{C}$, (i) $\bm{w}_{i}$ are conditionally independent given $X$; (ii) the sub-Gaussian norm of $\bm{w}_{i}$ is bounded by $K$ for all $i$. Let $\bm{W}:=[\bm{w}_{1},\bm{w}_{2},\dots,\bm{w}_{N}]{}^{\prime}$. Then for an $0<\epsilon<1$ we have

[TABLE]