# Implicit Bias of Projected Subgradient Method Gives Provable Robust Recovery
of Subspaces of Unknown Codimension

Paris Giampouras, Benjamin D. Haeffele & René Vidal
Mathematical Institute for Data Science
Johns Hopkins University
Baltimore, MD , USA
Correspondence to<EMAIL_ADDRESS>
###### Abstract

Robust subspace recovery (RSR) is a fundamental problem in robust
representation learning. Here we focus on a recently proposed RSR method
termed Dual Principal Component Pursuit (DPCP) approach, which aims to recover
a basis of the orthogonal complement of the subspace and is amenable to
handling subspaces of high relative dimension. Prior work has shown that DPCP
can provably recover the correct subspace in the presence of outliers, as long
as the true dimension of the subspace is known. We show that DPCP can provably
solve RSR problems in the unknown subspace dimension regime, as long as
orthogonality constraints -adopted in previous DPCP formulations- are relaxed
and random initialization is used instead of spectral one. Namely, we propose
a very simple algorithm based on running multiple instances of a projected
sub-gradient descent method (PSGM), with each problem instance seeking to find
one vector in the null space of the subspace. We theoretically prove that
under mild conditions this approach will succeed with high probability. In
particular, we show that 1) all of the problem instances will converge to a
vector in the nullspace of the subspace and 2) the ensemble of problem
instance solutions will be sufficiently diverse to fully span the nullspace of
the subspace thus also revealing its true unknown codimension. We provide
empirical results that corroborate our theoretical results and showcase the
remarkable implicit rank regularization behavior of PSGM algorithm that allows
us to perform RSR without being aware of the subspace dimension.

##  1 Introduction

Robust subspace recovery (RSR) refers to methods designed to identify an
underlying linear subspace (with dimension less than the ambient data
dimension) in a dataset which is potentially corrupted with outliers (i.e.,
points that do not lie in the linear subspace). Many methods for RSR have been
proposed in the literature over the past several years Xu et al. (2012); You
et al. (2017a); Lerman & Maunu (2018). Formulations based on convex
relaxations and decompositions of the data matrices into a low-rank matrices
plus sparse corruptions – either uniformly at random as in Robust PCA (RPCA)
Candès et al. (2011) or column-sparse corruptions as in Xu et al. (2012);
McCoy & Tropp (2011)111Note that RPCA is closely related but distinct from
RSR. In RPCA the corruptions are assumed to be entry-wise in a data matrix,
while in RSR the corruptions apply to an entire datapoint (i.e., column-wise).
– can, in certain situations, be shown to provably recover the true subspace
when the dimension is unknown. However, these theoretical guarantees often
require the dimension of the subspace, $d$, to be significantly less than the
ambient dimension of the data, $D$, and these methods are not suitable to the
more challenging high relative subspace dimension regime (i.e., when
$\frac{d}{D}\approx 1$).

The Dual Principal Component Pursuit approach for RSR. Recently, progress has
been made towards solving the RSR problem in the high relative dimension
regime by a formulation termed Dual Principal Component Pursuit (DPCP), which
is provably robust in recovering subspaces of high relative dimension,
Tsakiris & Vidal (2018). As implied by its name, DPCP follows a dual
perspective of RSR aiming to recover a basis for the orthogonal complement of
the inliers’ subspace. However, a key limitations of DPCP is that it requires
a priori knowledge of the true subspace dimension.

DPCP for $c=1$. Let $\tilde{{\bm{X}}}\in\mathbb{R}^{D\times(N+M)}$ denote the
data matrix defined as $\tilde{{\bm{X}}}=[{\bm{X}}\;\;{\bm{O}}]\bm{\Gamma}$
where ${\bm{X}}\in\mathbb{R}^{D\times N}$ is a matrix containing $N$ inliers
as its columns, ${\bm{O}}\in\mathbb{R}^{D\times M}$ is a matrix containing $M$
outliers, and $\bm{\Gamma}$ is an unknown permutation matrix. DPCP was first
formulated in Tsakiris & Vidal (2018) for handling subspaces of codimension
$c=D-d$ equal to 1 (i.e., the subspace is a hyperplane with dimension $d=D-1$)
formulated as the following optimization problem:

$\min_{{\bm{b}}\in\mathbf{R}^{D}}~{}~{}\|\tilde{{\bm{X}}}^{\top}{\bm{b}}\|_{1}\
\ \mathrm{s.t.}\ \ \|{\bm{b}}\|_{2}=1.$ (1)

Note that problem (1) is nonconvex due to the spherical constraint imposed on
vector ${\bm{b}}\in\mathbb{S}^{D-1}$ which is a normal vector of the $D-1$
dimensional hyperplane. In Tsakiris & Vidal (2018) the authors showed that the
global minimizer of (1) will give a normal vector of the underlying true
hyperplane as long as both inliers and outliers are well-distributed or the
ratio between the number of inliers and number of outliers is sufficiently
small. Following a probabilistic point of view, the authors in Zhu et al.
(2018) presented an improved theoretical analysis of DPCP giving further
insights on the remarkable robustness of DPCP in recovering the true
underlying subspaces even in datasets heavily corrupted by outliers. Moreover,
the authors introduced a projected subgradient method which converges to a
normal vector of the true subspace at a linear rate.

Recursive DPCP for known $c>1$. The authors of Zhu et al. (2018) also proposed
an extension to DPCP which allows for subspaces with codimension $c>1$ via a
projected subgradient algorithm which attempted to learn $c$ normal vectors to
the subspace in a recursive manner. Specifically, after convergence to a
normal vector, the projected subgradient algorithm is initialized with a
vector orthogonal to the previously estimated normal vector. However, for that
approach to be successful the knowledge of the true subspace codimension $c$
becomes critical. Specifically, if an underestimate of the true codimension
$c$ is assumed the recovered basis for the null space, $\hat{{\bm{B}}}$, will
fail to span the whole null space, $\mathcal{S}_{\perp}$. On the other hand,
an overestimate of $c$ will lead to columns of $\hat{{\bm{B}}}$ corresponding
to vectors that lie in $\mathcal{S}$.

Orthogonal DPCP for known $c$. In Zhu et al. (2019), an alternative to (1) was
proposed which attempts to solve for $c$ normal vectors to the subspace at
once via the formulation

$\min_{{\bm{B}}\in\mathbb{R}^{D\times
c}}\|\tilde{{\bm{X}}}^{\top}{\bm{B}}\|_{1,2}~{}~{}~{}\ \ \mathrm{s.t.}\ \
~{}~{}{\bm{B}}^{\top}{\bm{B}}=\mathbf{I}.$ (2)

The authors also propose an optimization algorithm based on the projected
Riemannian subgradient method (RSGM), which builds on similar ideas as the
projected subgradient method of Zhu et al. (2018) and enjoys a linear converge
rate when the step size is selected based on a geometrically diminishing rule.
In Ding et al. (2021) a theoretical analysis is provided on the geometric
properties of 2 showing the merits of DPCP in handling datasets a) highly
contaminated by outliers (in the order of $M=\mathcal{O}(N^{2})$) and b)
subspaces of high relative dimension. However again, it is critical to note
that a key shortcoming of this approach is the fact that because all
minimizers of problem (2) will be orthogonal matrices, a prerequisite for
recovering the correct orthogonal complement of the inliers subspace is the a
priori knowledge of the true codimension $c$ (see Fig. 1).

Contributions. In this work, we address this key limitation by proposing a
framework that allows us to perform robust subspace recovery in high relative
subspace dimension regime without requiring a priori knowledge of the true
subspace dimension. In particular, our proposed approach will be based on the
simple approach of solving multiple, parallel instances of the DPCP
formulation for solving for a single normal vector to the subspace,

$\min_{{\bm{B}}\in\mathbb{R}^{D\times
c^{\prime}}}~{}\sum^{c^{\prime}}_{i=1}\|\tilde{{\bm{X}}}^{\top}{\bm{b}}_{i}\|_{1}~{}\
\ \mathrm{s.t.}\ \
~{}~{}\|{\bm{b}}_{i}\|_{2}=1,~{}~{}~{}i=1,2,\dots,c^{\prime}$ (3)

where $c^{\prime}$ is assumed to be a safe overestimate of $c$ i.e.,
$c^{\prime}\geq c$. Contrary to (2), the objective function in (3) decouples
over the columns ${\bm{b}}_{i}$ of matrix
${\bm{B}}=[{\bm{b}}_{1}~{}{\bm{b}}_{2}~{}\cdots{\bm{b}}_{c^{\prime}}]$ and
thus can be can be solved in a parallel manner by independently applying a
projected subgradient algorithm (referred to as PSGM) from $c^{\prime}$
different random initialization. Moreover, we observe that with random
initialization we can get vectors sufficiently spread on the sphere that lead
PSGM (initialized with those vectors) to return normal vectors of
$\mathcal{S}$. These are all linearly independent when $c^{\prime}\leq c$ and
thus can span $\mathcal{S}_{\perp}$ when $c^{\prime}=c$. If $c^{\prime}>c$
then PSGM will return $c^{\prime}-c$ redundant vectors that will still lie in
$\mathcal{S}_{\perp}$ yet they will be linearly dependent (see Figure 1). That
being said, we show that this simple strategy permits us to robustly recover
the true subspace even without knowledge of the true codimension $c$.

As is detailed in Sections 3 and 4, this remarkable behavior of PSGM
originates from the implicit bias that is induced in the optimization problem
due to a) the relaxation of orthogonality constraints in (3) and b) the random
initialization scheme that is adopted. Our specific contributions are as
follows:

  1. 1.

First, we focus on a continuous version of (3) i.e., in the case that inliers
and outliers are distributed under continuous measures which induces a benign
landscape in DPCP and hence is easier to be analyzed. We prove that DPCP
problem formulated as in (3) can be solved via a projected subgradient
algorithm that implicitly biases solutions towards low-rank matrices
$\hat{{\bm{B}}}\in\mathbb{R}^{D\times c^{\prime}}$ whose columns are the
projections of the randomly initialized columns of ${\bm{B}}^{0}$ onto
$\mathcal{S}_{\perp}$. As a result, $\hat{{\bm{B}}}$ almost surely spans
$\mathcal{S}_{\perp}$ as long as it is randomly initialized with
$c^{\prime}\geq c$.

  2. 2.

Second, we analyze the discrete version which is more challenging, yet of more
practical interest, showing that iterates of DPCP-PSGM converge to a scaled
and perturbed version of the initial matrix ${\bm{B}}^{0}$. This compelling
feature of DPCP-PSGM allows to derive a sufficient condition and a
probabilistic bound guaranteeing when the matrix
$\hat{{\bm{B}}}\in\mathbb{R}^{D\times c^{\prime}}$ spans
$\mathcal{S}_{\perp}$.

  3. 3.

We provide empirical results both on simulated and a real datasets,
corroborating our theory and showing the robustness of our approach even
without knowledge of the true subspace codimension.

|
---|---
Figure 1: Graphical illustration of the recovered normal vectors of
$\mathcal{S}$ by (left) the proposed DPCP-PSGM approach and (right) methods
that use spectral initialization and impose orthogonality constraints .
Initial vectors ${\bm{b}}^{0}_{1},{\bm{b}}^{0}_{2},{\bm{b}}^{0}_{3}$ are
randomly initialized and are non-orthogonal in (left) and spectrally
initialized (hence orthogonal) in (right). Note that in (left)
$\operatorname*{rank}(\hat{{\bm{B}}}^{\ast})$ (where
$\hat{{\bm{B}}}^{\ast}=[{\bm{b}}^{\ast}_{1},{\bm{b}}^{\ast}_{2},{\bm{b}}^{\ast}_{3}]$)
equals to the true codimension $c=2$ of $\mathcal{S}$ and
$\operatorname*{span}({\bm{B}}^{\ast})\equiv\mathcal{S}_{\perp}$ while in
(right) ${\bm{B}}^{\ast}$ is orthogonal hence
$\operatorname*{rank}({\bm{B}}^{\ast})=3$ with
${\bm{b}}^{\ast}_{2}\in\mathcal{S}$.

##  2 Related work

Subspace Recovery. Learning underlying low-dimensional subspace
representations of data has been a central topic of interest in machine
learning research. Principal Component Analysis (PCA) has been the most
celebrated method of this kind and is based on the minimization of the
perpendicular distances of the data points from the estimates linear subspace,
Jolliffe & Cadima (2016). Albeit, it is originally formulated as nonconvex
optimization problem, PCA can be easily solved in closed form using a singular
value decomposition (SVD) operation (see e.g. Vidal et al. (2016)). Despite
its great success, PCA is prone to failure when handling datasets that contain
outliers i.e., data points whose deviation from the inliers’ subspace is
“large” in the $\ell_{2}$ norm sense.

Robust Subspace Recovery (RSR). To remedy this weakness of PCA robust subspace
recovery (RSR) methods attempt to identify the outliers in the dataset an
recover the true underlying low-dimensional subspace of the inliers Lerman &
Maunu (2018); Maunu et al. (2019). A classical approach to this problem is
RANSAC Fischler & Bolles (1981), which is given a time budget and randomly
chooses per iteration $d$ points and then fits a $d$-dimensional subspace to
those points and checks how the proposed subspace fits the remaining data
points. RANSAC then outputs the subspace that agrees with the largest number
of points. However, RANSAC’s reliance on randomly sampling points to propose
subspaces can be highly inefficient when the number of outliers is high (as
well as the fact that RANSAC also needs knowledge of the true subspace
dimension, $d$). The need to tackle inherent shortcomings of RANSAC pertaining
to computational complexity issues inspired alternative convex formulations of
RSR, Xu et al. (2012); You et al. (2017b); Rahmani & Atia (2017); Zhang &
Lerman (2014). In Xu et al. (2012) the authors decompose the data matrix as a
sum of a low-rank and a column-sparse component. However, theoretical
guarantees obtained for convex formulations only hold for subspaces of
relatively low-dimensional subspaces i.e., for $d\ll D$ where $d$ and $D$
denote the subspace and the ambient dimension, respectively. To the best of
our knowledge, existing RSR algorithms rely heavily on one of two key
assumptions. 1) The subspace is very low-dimensional relative to the ambient
dimension $(d\ll D)$ or 2) The subspace dimension is a priori known.
Undoubtedly, the second hypothesis is rather strong in real-world
applications, and many applications also do not satisfy the first assumption.
Moreover, heuristic strategies for selecting the dimension of the subspace are
hard to be applied in the RSR setting since they incur computationally
prohibitive procedures, Lerman & Maunu (2018).
Relation to Orthogonal Dictionary Learning (ODL). Note that objective
functions in the form of (3) show up beyond RSR problems i.e., in orthogonal
dictionary learning (ODL), sparse blind deconvolution, etc., Qu et al. (2020).
Specifically, based on a similar formulation the authors in Bai et al. (2019)
proved that $c^{\prime}=\mathcal{O}(c\log c)$ independent random initial
vectors suffice in order to recover with high probability a dictionary of size
$D\times c$ with high accuracy. In this paper we aim to recover a basis of the
orthogonal complement of a subspace of unknown dimension instead of accurately
estimating a dictionary hence our goal differs from that in Bai et al. (2019).
Implicit bias in Robust Recovery Problems. The notions of implicit bias and
implicit regularization have been used interchangeably in the nonconvex
optimization literature for describing the tendency of optimization algorithms
to converge to global minima of minimal complexity with favorable
generalization properties in overparameterized models, Gunasekar et al.
(2018). In the context of robust recovery, the authors in You et al. (2020)
showed that Robust PCA can be suitably re-parametrized in such a way to favor
low-rank and sparse solutions without using any explicit regularization. In
this work, we use the term implicit bias for describing the convergence of
DPCP-PSGM to low-rank solutions, which are not necessarily global minimizers,
that span the orthogonal complement of the subspace when a) orthogonality
constraints in DPCP formulation are relaxed b) DPCP is overparameterized i.e.,
$c^{\prime}\geq c$ and c) PSGM randomly initialized. To the best of our
knowledge our work is the first that explores implicit bias in projected
subgradient algorithms that minimize objective functions with
sphere/orthogonality constraints.

##  3 Dual Principal Component Pursuit and the Projected Subgradient Method

We re-write the DPCP formulation given in (1) as

$\min_{{\bm{B}}\in\mathbb{R}^{D\times
c^{\prime}}}~{}\|\tilde{{\bm{X}}}^{\top}{\bm{b}}\|_{1}=\|{\bm{X}}^{\top}{\bm{b}}\|_{1}+\|{\bm{O}}^{\top}{\bm{b}}\|_{1}~{}\
\ \mathrm{s.t.}\ \ ~{}~{}\|{\bm{b}}\|_{2}=1$ (4)

In Zhu et al. (2018), the authors proposed a projected subgradient descent
algorithm for addressing (4) that consists of a subgradient step followed by a
projection onto the sphere i.e.,

$\begin{split}{\bm{b}}^{k+1}=\hat{{\bm{b}}}^{k}-\mu^{k}\left({\bm{X}}\mathrm{Sgn}({\bm{X}}^{\top}{\bm{b}}^{k})+{\bm{O}}\mathrm{Sgn}({\bm{O}}^{\top}{\bm{b}}^{k})\right)\
\ \text{and}\ \
\hat{{\bm{b}}}^{k+1}=\mathcal{P}_{\mathbb{S}^{D-1}}({\bm{b}}^{k+1}),\end{split}$
(5)

where $\mu^{k}$ is the -adaptively updated per iteration- step size and
$\hat{{\bm{b}}}^{k}$ is the unit $\ell_{2}$ norm vector corresponding to the
$k$th iteration.

The convergence properties of the projected subgradient algorithm described
above depend on specific quantities denoted as $c_{{\bm{X}},\min}$ and
$c_{{\bm{X}},\max}$ that reflect the geometry of the problem and are defined
as
$c_{{\bm{X}},\min}=\frac{1}{N}\min_{{\bm{b}}\in\mathbb{S}^{d-1}\cap\mathcal{S}}\|{\bm{X}}^{\top}{\bm{b}}\|_{1}$
and
$c_{{\bm{X}},\max}=\frac{1}{N}\max_{{\bm{b}}\in\mathbb{S}^{D-1}\cap\mathcal{S}}\|{\bm{X}}^{\top}{\bm{b}}\|_{1}$.
Note that the more well distributed the inliers are in the subspace
$\mathcal{S}$ the higher the value of the quantity $c_{{\bm{X}},\min}$ (called
as permeance statistic which first appeared in Lerman et al. (2015)) as it
becomes harder to find a vector ${\bm{b}}$ in the subspace $\mathcal{S}$ that
is orthogonal to the inliers. Moreover, $c_{{\bm{X}},\min}$ and
$c_{{\bm{X}},\max}$ converge to the same value as $N\rightarrow\infty$
provided the inliers are uniformly distributed in the subspace i.e.,
$c_{{\bm{X}},\min}\rightarrow c_{d},c_{{\bm{X}},\max}\rightarrow c_{d}$, where
$c_{d}$ is given as the average height of the unit hemisphere on
$\mathbb{R}^{d}$,

$c_{d}:=\frac{(d-2)!!}{(d-1)!!}\Bigg{\\{}\begin{matrix}\frac{2}{\pi},~{}\text{if}~{}~{}d~{}~{}\text{is
~{}even},\\\\[5.69046pt] 1,~{}\text{if}~{}d~{}~{}\text{is
odd}~{}\end{matrix}~{}~{}~{}\text{where
}k!!=\Bigg{\\{}\begin{matrix}k(k-2)(k-4)\cdots 4\cdot 2,~{}k~{}~{}\text{is
~{}even},\\\\[5.69046pt] k(k-2)(k-4)\cdots 3\cdot 1,~{}k~{}~{}\text{is
odd}~{}\end{matrix}$ (6)

Similarly to $c_{{\bm{X}},\min},c_{{\bm{X}},\max}$, we will also be interested
in quantities $c_{{\bm{O}},\min},c_{{\bm{O}},\max}$ which indicate how well-
distributed the outliers are in the ambient space. These quantities are
defined as
$c_{{\bm{O}},\min}=\min_{{\bm{b}}\in\mathbb{S}^{D-1}}\frac{1}{M}\|{\bm{O}}^{\top}{\bm{b}}\|_{1}$
and
$c_{{\bm{O}},\max}=\max_{{\bm{b}}\in\mathbb{S}^{D-1}}\frac{1}{M}\|{\bm{O}}^{\top}{\bm{b}}\|_{1}$.
$c_{{\bm{O}},\max}$ can be viewed as the dual permeance statistic and is
bounded away from small values while its difference from $c_{{\bm{O}},\min}$
tends to zero as $M\rightarrow\infty$. Further, if the outliers are uniformly
distributed on the sphere, then $c_{{\bm{O}},\max}\rightarrow c_{D}$ and
$c_{{\bm{O}},\min}\rightarrow c_{D}$ where $c_{D}$ is defined as in (6), (Zhu
et al., 2018).

Finally, we also define the quantities
$\eta_{{\bm{O}}}=\frac{1}{M}\max_{{\bm{b}}\in\mathbb{S}^{D-1}}\|({\bm{I}}-{\bm{b}}{\bm{b}}^{\top}){\bm{O}}\mathrm{Sgn}({\bm{O}}^{\top}{\bm{b}})\|_{2}$
and
$\eta_{{\bm{X}}}=\frac{1}{M}\max_{{\bm{b}}\in\mathbb{S}^{D-1}}\|(\mathcal{P}_{\mathcal{S}}-{\bm{b}}{\bm{b}}^{\top}){\bm{X}}\mathrm{Sgn}({\bm{X}}^{\top}{\bm{b}})\|_{2}$.
As $M\rightarrow\infty$ and assuming outliers in ${\bm{O}}$ are well-
distributed we get ${\bm{O}}\mathrm{Sgn}({\bm{O}}^{\top}{\bm{b}})\rightarrow
c_{D}{\bm{b}}$ thus $\eta_{{\bm{O}}}\rightarrow 0$ (Tsakiris & Vidal, 2018).
Likewise, $\eta_{{\bm{X}}}\rightarrow 0$ as $N\rightarrow\infty$ provided that
inliers are uniformly distributed in the $d$-dimensional subspace. The
following theorem (see full version in Appendix) provides convergence
guarantees of the projected subgradient method that was proposed in Zhu et al.
(2018) for addressing problem (1).

###### Theorem 1

(Informal Theorem 3 of Zhu et al. (2018)) Let $\\{\hat{{\bm{b}}}_{k}\\}$ the
sequence generated by the projected subgradient algorithm in Zhu et al.
(2018), with initialization $\hat{{\bm{b}}}_{0}$ such that

$\theta_{0}<\mathrm{arctan}\left(\frac{Nc_{{\bm{X}},\min}}{N\eta_{{\bm{X}}}+M\eta_{{\bm{O}}}}\right)\
\ \textrm{and}\ \ Nc_{{\bm{X}},\min}\geq N\eta_{{\bm{X}}}+M\eta_{{\bm{O}}}$
(7)

where $\theta_{0}$ denotes the principal angle of ${\bm{b}}^{0}$ from
$\mathcal{S}_{\perp}$. If the step size $\mu^{k}$ is updated according to a
piecewise geometrically diminishing rule given as

$\mu^{k}=\left\\{\begin{array}[]{l l}\mu^{0},&k<K_{0}\\\
\mu^{0}\beta^{\lfloor(k-K_{0})/K_{\ast}\rfloor+1},&k\geq
K_{0}\end{array}\right.$ (8)

where $\beta<1$, $\lfloor\cdot\rfloor$ is the floor function, then the
iterates ${\bm{b}}^{k}$ converge to a normal vector of $\mathcal{S}$.

##  4 Dual Principal Component Pursuit in subspaces of unknown codimension

Current theoretical results provide guarantees for recovering the true inlier
subspace, when the proposed algorithms know a priori of the subspace
codimension $c$, which is a rather strong requirement and is far from being
true in real word applications. Here we describe our proposed approach, which
consists of removing the orthogonality constraint on ${\bm{B}}$, along with a
theoretical analysis that gives guarantees of recovering the true underlying
subspace even when the true codimension $c$ is unknown. First we analyze a
continuous version of DPCP, which arises when the number of inliers and
outliers are distributed according to continuous measures and their number
tends to $\infty$. The continuous DPCP incurs an optimization problem with a
benign landscape that allows us to better illustrate the favorable properties
of DPCP-PSGM when it comes to the convergence of its iterates. Then we extend
the results to the discrete case that deals with a finite number of inliers
and outliers yielding a more challenging optimization landscape.

###  4.1 PSGM’s iterates convergence in the continuous version of DPCP

The following lemma provides the continuous version of the discrete objective
function given in (3).

###### Lemma 2

In the continuous case, the discrete DPCP problem given in (3) is reformulated
as,

$\begin{split}\min_{{\bm{B}}\in\mathbb{R}^{D\times
c^{\prime}}}\sum^{c^{\prime}}_{i=1}\left(p\mathbb{E}_{\bm{\mu}_{\mathbb{S}^{D-1}}}[f_{\mathbf{b}_{i}}]+(1-p)\mathbb{E}_{\bm{\mu}_{\mathbb{S}^{D-1}\cap\bm{S}}}[f_{{\bm{b}}_{i}}]\right)=\sum^{c^{\prime}}_{i=1}\|{\bm{b}}_{i}\|_{2}\left(pc_{D}+(1-p)c_{d}\mathrm{cos}(\phi_{i})\right)\\\
\ \ \mathrm{s.t.}\ \
\|{\bm{b}}_{i}\|_{2}=1,~{}~{}i=1,2,\dots,c^{\prime}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(9)

where $f_{{\bm{b}}}:\mathbb{S}^{D-1}\rightarrow\mathbb{R}$,
$f_{{\bm{b}}}(\mathbf{z})=|\mathbf{z}^{\top}{\bm{b}}|$, $\phi_{i}$ is the
principal angle of ${\bm{b}}_{i}$ from the inliers subspace $\mathcal{S}$ and
$p$ is the probability of occurrence of an outlier.

Note that $\bm{\mu}_{\mathbb{S}^{D-1}},\bm{\mu}_{\mathbb{S}^{D-1}\cap\bm{S}}$
are the continuous measures associated with the outliers and inliers,
respectively. Evidently, (9) attains its global minimum for vectors
${\bm{b}}_{i}$s that are orthogonal to the inliers’ subspace. Based on (3) and
due to Lemma 2, we can now minimize the objective function of the “continuous
version” of DPCP by employing a projected subgradient methods (PSGM) that
performs the following steps per iteration 222Note that
$\partial\|{\bm{b}}\|_{2}=\frac{{\bm{b}}}{\|{\bm{b}}\|_{2}}$ for
${\bm{b}}\neq\mathbf{0}$ and
$\|{\bm{b}}_{i}\|_{2}\mathrm{cos}(\phi_{i})={\bm{b}}^{\top}_{i}\hat{\mathbf{s}}_{i}$
where
$\hat{\mathbf{s}}_{i}=\frac{\mathcal{P}_{\mathcal{S}}({\bm{b}}_{i})}{\|\mathcal{P}_{\mathcal{S}}({\bm{b}}_{i})\|_{2}}$.

$\begin{split}{\bm{b}}^{k+1}_{i}=\hat{{\bm{b}}}^{k}_{i}-\mu^{k}_{i}(pc_{D}\hat{{\bm{b}}}^{k}_{i}+(1-p)c_{d}\hat{\mathbf{s}}^{k}_{i})\
\ \text{and}\ \
\hat{{\bm{b}}}^{k+1}_{i}=\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{k+1}_{i})~{},i=1,2,\dots,c^{\prime}\end{split}$
(10)

###### Lemma 3

A projected subgradient algorithm consisting of the steps described in (10)
using a piecewise geometrically diminishing step size rule (see (8) in Theorem
1) will almost surely asymptotically converge to a matrix
$\hat{{{\bm{B}}}}^{\ast}\in\mathbb{R}^{D\times c^{\prime}}$ whose columns
$\hat{{\bm{b}}}^{\ast}_{i},~{}i=1,2,\dots,c^{\prime}$ will be normal vectors
of the inliers’ subspace when randomly initialized with vectors
${\bm{b}}^{0}_{i}\in\mathbb{S}^{D-1},~{}i=1,2,\dots,c^{\prime}$ uniformly
distributed over the sphere $\mathbb{S}^{D-1}$.

Lemma 3 allows us to claim that we can always recover $c^{\prime}\geq c$
normal vectors to the inliers’ subspace using a PSGM algorithm consisting of
steps given in (10). However, this does not tell the whole story yet, since
our ultimate objective is to recover a matrix $\hat{{\bm{B}}}$ that spans
$\mathcal{S}_{\perp}$. Thus, it remains to show that the rank of
$\hat{{\bm{B}}}$ is equal to the true and unknown codimension of the inliers’
subspace $c$. Next we prove that by initializing with a $\hat{{\bm{B}}}_{0}$
such that $\operatorname*{rank}(\hat{{\bm{B}}}_{0})=c^{\prime}$ (i.e.,
$\hat{{\bm{B}}}_{0}$ is initialized to be full-rank), we can guarantee that we
can solve the continuous version of DPCP using PSGM and converge to a
$\hat{{\bm{B}}}$ such that $\operatorname*{rank}(\hat{{\bm{B}}})=c$ thus
getting $\operatorname*{span}(\hat{{\bm{B}}})\equiv\mathcal{S}_{\perp}$ (along
with recovering the true subspace dimension). By projecting the PSGM iterates
given in (10) onto $\mathcal{S}_{\perp}$ we have,

$\begin{split}\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{k+1}_{i})=(1-\mu^{k}_{i}pc_{D})\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{k}_{i})\
\ \text{and}\ \
\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{k+1}_{i})=\mathcal{P}_{\mathcal{S}_{\perp}}(\mathcal{P}_{\mathbb{S}^{D-1}}({\bm{b}}^{k+1}_{i}))\end{split}$
(11)

We hence observe that the projections of successive iterates of PSGM are
scaled versions of the corresponding projections of the previous iterates. We
can now state Lemma 4.

###### Lemma 4

The PSGM iterates
$\hat{{\bm{b}}}^{k}_{i},~{}i=1,2,\dots,c^{\prime},~{}k=1,2,\dots$ given in
(10), when randomly initialized with $\hat{{\bm{b}}}^{0}_{i}$s,
$i=1,2,\dots,c^{\prime}$ that are independently drawn from a spherical
distribution with unit $\ell_{2}$ norm converge almost surely to $c^{\prime}$
normal vectors of the inliers subspace $\mathcal{S}$ denoted as
$\hat{{\bm{b}}}^{\ast}_{i},~{}i=1,2,\dots,c^{\prime}$ that are given by
$\hat{{\bm{b}}}^{\ast}_{i}=\frac{\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{i})}{\|\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{i})\|_{2}},~{}~{}~{}~{}~{}i=1,2,\dots,c^{\prime}$.

Lemma 4 shows that the initialization of PSGM plays a pivotal role since it
determines the direction of the recovered normal vectors
$\\{\hat{{\bm{b}}}^{\ast}_{i}\\}_{i=1}^{c^{\prime}}$. Lemmas 3 and 4 pave the
way for Theorem 5.

###### Theorem 5

Let $\hat{{\bm{B}}^{0}}\in\mathbb{R}^{D\times c^{\prime}}$ where
$c^{\prime}\geq c$ with $c$ denoting the true codimension of the inliers
subspace $\mathcal{S}$, consisting of unit $\ell_{2}$ norm column vectors
$\hat{{\bm{b}}}^{0}_{i}\in\mathbb{S}^{D-1},i=1,2,\dots,c^{\prime}$ that are
independently drawn from uniform distribution over the sphere
$\mathbb{S}^{D-1}$. A PSGM algorithm initialized with $\hat{{\bm{B}}^{0}}$
will almost surely converge to a matrix $\hat{{\bm{B}}}^{\ast}$ such that
$\operatorname*{span}(\hat{{\bm{B}}}^{\ast})\equiv\mathcal{S}_{\perp}$.

From Theorem 5 we observe that in the benign scenario where inliers and
outliers are distributed under continuous measures, we can recover the correct
orthogonal complement of the inlier’s subspace even when we are oblivious to
its true codimension. Remarkably, this is achieved by exploiting the implicit
bias induced by multiple random initializations of the PSGM algorithm for
solving the DPCP formulation given in (3), which is free of orthogonality
constraints.

###  4.2 PSGM’s iterates convergence in the discrete version of DPCP

From this analysis of the continuous version of DPCP we now extend to the the
discrete version, which is of more practical relevance for finite data, yet
also presents more challenges. To begin, we reformulate the DPCP objective as
follows

$\sum^{c^{\prime}}_{i=1}\|\tilde{{\bm{X}}}^{\top}{\bm{b}}_{i}\|_{1}=\sum^{c^{\prime}}_{i=1}\|{\bm{X}}^{\top}{\bm{b}}_{i}\|_{1}+\|{\bm{O}}^{\top}{\bm{b}}_{i}\|_{1}=M\sum^{c^{\prime}}_{i=1}\mathbf{b}^{\top}_{i}\bm{\mathit{o}}_{{\bm{b}}_{i}}+N\sum^{c^{\prime}}_{i=1}{\bm{b}}^{\top}_{i}\bm{\mathit{x}}_{{\bm{b}}_{i}}$
(12)

where $\bm{\mathit{x}}_{{\bm{b}}_{i}}$ and $o_{{\bm{b}}_{i}}$ are called as
average inliers and average outliers terms, defined as
$\bm{\mathit{x}}_{{\bm{b}}_{i}}=\frac{1}{N}\sum^{N}_{j=1}Sgn({\bm{b}}^{\top}_{i}\bm{\mathit{x}}_{j})\bm{\mathit{x}}_{j}$
and
$\bm{\mathit{o}}_{{\bm{b}}_{i}}=\frac{1}{M}\sum^{M}_{j=1}Sgn({\bm{b}}^{\top}_{i}\bm{\mathit{o}}_{j})\bm{\mathit{o}}_{j}$.

In Algorithm 1, we give the projected subgradient method (DPCP-PSGM) applied
on the DPCP problem given in (3).

Result:
$\hat{{\bm{B}}}=[\hat{{\bm{b}}}^{k}_{1},\hat{{\bm{b}}}^{k}_{2},\dots,\hat{{\bm{b}}}^{k}_{c^{\prime}}]$

Initialize: Randomly sample
$\hat{{\bm{b}}}^{1}_{0},\hat{{\bm{b}}}^{2}_{0},\dots,\hat{{\bm{b}}}^{c^{\prime}}_{0}$
from a uniform distribution on $\mathbb{S}^{D-1}$ ;

for  _$k=1,2,\dots$_ do

for  _$i=1,2,\dots,c^{\prime}$_ do

Update the step-size according to a specific rule;

${\bm{b}}^{k+1}_{i}=\hat{{\bm{b}}}^{k}_{i}-\mu^{k}_{i}(M\bm{\mathit{o}}^{k}_{\hat{{\bm{b}}}_{i}}+N\bm{\mathit{x}}^{k}_{\hat{{\bm{b}}}_{i}})$;

$\hat{{\bm{b}}}^{k+1}_{i}=\mathcal{P}_{\mathbb{S}^{D-1}}({\bm{b}}^{k+1}_{i})$;

end for

end for

Algorithm 1 DPCP-PSGM algorithm for solving (3)

Note that the average outliers and inliers terms are discrete versions of the
corresponding continuous average terms $c_{D}{\bm{b}}_{i}$ and
$c_{d}\hat{\mathbf{s}}_{i}$ where
$\hat{\mathbf{s}}_{i}=\mathcal{P}_{\mathcal{S}}({\bm{b}}_{i})$, respectively,
Tsakiris & Vidal (2018). We now express the sub-gradient step of Algorithm 1
as

${\bm{b}}^{k+1}_{i}=\hat{{\bm{b}}}^{k}_{i}-\mu^{k}_{i}\left(M(c_{D}\hat{{\bm{b}}}^{k}_{i}+\bm{\mathit{e}}_{{\bm{O}}}^{i,k})+N(c_{D}\hat{\mathbf{s}}^{k}_{i}+\bm{\mathit{e}}_{{\bm{X}}}^{i,k})\right)$
(13)

where the quantities
$\bm{\mathit{e}}_{{\bm{O}}}^{i,k}=\bm{\mathit{o}}_{\hat{{\bm{b}}}^{k}_{i}}-c_{D}\hat{{\bm{b}}}^{k}_{i}-$
and
$\bm{\mathit{e}}_{{\bm{X}}}^{i,k}=\bm{\mathit{x}}_{\hat{{\bm{b}}}^{k}_{i}}-c_{D}\hat{\mathbf{s}}^{k}_{i}$
account for the error between the continuous and discrete versions of the
average outliers and the average inliers terms, respectively. Following a
similar path as in the continuous case we next project the iterates of (13)
onto $\mathcal{S}_{\perp}$,

$\begin{split}\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{k+1}_{i})&=\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{k}_{i})-\mu^{k}_{i}\left(M\mathcal{P}_{\mathcal{S}_{\perp}}(c_{D}\hat{\mathbf{b}}^{k}_{i}+\bm{\mathit{e}}_{{\bm{O}}}^{i,k})+N\cancelto{0}{\mathcal{P}_{\mathcal{S}_{\perp}}(c_{D}\hat{\mathbf{s}}^{k}_{i}+\bm{\mathit{e}}_{{\bm{X}}}^{i,k})}\right)\\\
&=(1-\mu^{k}_{i}Mc_{D})\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{k}_{i})-\mu^{k}_{i}M\mathcal{P}_{\mathcal{S}_{\perp}}(\bm{\mathit{e}}_{{\bm{O}}}^{i,k})\end{split}$
(14)

Remark. Eq. (14) reveals that DPCP-PSGM, applied on the discrete problem,
gives rise to updates whose projections to $\mathcal{S}_{\perp}$ are scaled
and perturbed versions of the previous estimates. The magnitude of
perturbation depends on the discrepancy between the continuous and the
discrete problem.

Recall that
$\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{k}_{i})=\frac{\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{k}_{i})}{\|{\bm{b}}^{k}_{i}\|_{2}}$,
and we can rewrite the update of the 2nd iteration of DPCP-PSGM,

$\begin{split}&\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{2}_{i})=\frac{(1-\mu^{1}_{i}Mc_{D})}{\|{\bm{b}}^{1}_{i}\|_{2}}\left((1-\mu^{0}_{i}Mc_{D})\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{i})-\mu^{0}_{i}M\mathcal{P}_{\mathcal{S}_{\perp}}(\bm{\mathit{e}}_{{\bm{O}}}^{i,0})\right)-\mu^{1}_{i}M\mathcal{P}_{\mathcal{S}_{\perp}}(\bm{\mathit{e}}_{{\bm{O}}}^{i,1})\\\
&=\frac{(1-\mu^{1}_{i}Mc_{D})(1-\mu^{0}_{i}Mc_{D})}{\|{\bm{b}}^{1}_{i}\|_{2}\|{\bm{b}}^{0}_{i}\|_{2}}\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{\mathbf{b}}^{0}_{i})-\frac{(1-\mu^{1}_{i}Mc_{D})}{\|{\bm{b}}^{1}_{i}\|_{2}}\mu^{0}_{i}M\mathcal{P}_{\mathcal{S}_{\perp}}(\bm{\mathit{e}}_{{\bm{O}}}^{i,0})-\mu^{1}_{i}M\mathcal{P}_{\mathcal{S}_{\perp}}(\bm{\mathit{e}}_{{\bm{O}}}^{i,1})\end{split}$
(15)

where we have assumed that $\|{\bm{b}}^{0}_{i}\|_{2}=1$. By repeatedly
applying the same steps, we can reach to the following recursive expression
for $\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{K}_{i})$,

$\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{K}_{i})=\left(\prod^{K-1}_{k=0}\frac{(1-\mu^{k}_{i}Mc_{D})}{\|{\bm{b}}^{k}_{i}\|_{2}}\right)\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{0}_{i})-\sum^{K-1}_{k=0}\left(\prod^{K-1}_{j=k+1}\frac{(1-\mu^{j}_{i}Mc_{D})}{\|{\bm{b}}^{j}_{i}\|_{2}}\right)\mu^{k}_{i}M\mathcal{P}_{\mathcal{S}_{\perp}}(\bm{\mathit{e}}_{{\bm{O}}}^{i,k})$
(16)

where for $j>K-1$ we set
$\prod^{K-1}_{j=k+1}\frac{(1-\mu^{j}_{i}Mc_{D})}{\|{\bm{b}}^{j}_{i}\|_{2}}=1$.

By dividing (16) with
$\prod^{K-1}_{k=0}\frac{(1-\mu^{k}_{i}Mc_{D})}{\|{\bm{b}}^{k}_{i}\|_{2}}$ and
by projecting onto the sphere $\mathbb{S}^{D-1}$ we get

$\mathcal{P}_{\mathbb{S}^{D-1}}(\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{K}_{i}))=\mathcal{P}_{\mathbb{S}^{D-1}}(\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{0}_{i})-\mathcal{P}_{\mathcal{S}_{\perp}}(\bm{\delta}^{K}_{i}))$
(17)

where $\bm{\delta}_{i}$ is defined as
$\bm{\delta}^{K}_{i}=\sum^{K-1}_{k=0}\left(\prod^{k}_{j=0}\frac{\|{\bm{b}}^{j}_{i}\|_{2}}{(1-\mu^{j}_{i}Mc_{D})}\right)\mu^{k}_{i}M\bm{\mathit{e}}_{{\bm{O}}}^{i,k}$.

Assumption 1. We assume that the principal angles $\theta^{i}_{0}$ for all
${\bm{b}}^{0}_{i}$s satisfy the inequality
$\theta^{i}_{0}<\mathrm{arctan}\left(\frac{Nc_{{\bm{X}},\min}}{N\eta_{{\bm{X}}}+M\eta_{{\bm{O}}}}\right)~{}~{}\forall
i,i=1,2,\dots,c^{\prime}$.

Assumption 1 essentially assumes that the sufficient condition given in eq.
(7) required by PSGM algorithm for converging to a normal vector is satisfied
which is the same condition for success in Zhu et al. (2018). Under Assumption
1 we can invoke the convergence properties of PSGM given in Theorem 1 and get
as
$K\rightarrow\infty,\mathcal{P}_{\mathbb{S}^{D-1}}(\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{K}))\rightarrow\hat{{\bm{b}}}^{\ast}\in\mathcal{S}_{\perp}\cap\mathbb{S}^{D-1}$.
That being said, we denote
$\hat{{\bm{b}}}^{\ast}_{i}=\mathcal{P}_{\mathbb{S}^{D-1}}\left(\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{0}_{i})-\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{\bm{\delta}_{i}})\right)$,
where
$\hat{\bm{\delta}_{i}}=\lim_{K\rightarrow\infty}\bm{\delta}^{K}_{i}$333For the
sake of brevity we assume that the step size has been selected such that
existence of the limit is guaranteed. We refer the reader to the proof Lemma 8
for further details.. Following the same steps as in Section 4.1 and by
defining matrices
${\bm{B}}^{\ast}=[{\bm{b}}^{\ast}_{1},{\bm{b}}^{\ast}_{2},\dots,{\bm{b}}^{\ast}_{c^{\prime}}],{\bm{B}}^{0}=[{\bm{b}}^{0}_{1},{\bm{b}}^{0}_{2},\dots,{\bm{b}}^{0}_{c^{\prime}}],\hat{\bm{\Delta}}=[\hat{\bm{\delta}}_{1},\hat{\bm{\delta}}_{2},\dots,\hat{\bm{\delta}}_{c^{\prime}}]$
we can express the matrix ${\bm{B}}^{\ast}$ as
$\hat{{\bm{B}}}^{\ast}=\mathcal{P}_{\mathbb{S}^{D-1}}\left(\mathcal{P}_{\mathcal{S}_{\perp}}\left({\bm{B}}^{0}-\hat{\bm{\Delta}}\right)\right)$
where ${\bm{B}}^{\ast}$ will now consist of normal vectors of the inliers’
subspace. In order to guarantee that
$\operatorname*{span}({{\bm{B}}^{\ast}})\equiv\mathcal{S}_{\perp}$ it thus
suffices to ensure that $\operatorname*{rank}\left({\bm{B}}^{\ast}\right)=c$.
Here we show that a sufficient condition for this to hold is that the matrix
${\bm{A}}={\bm{B}}^{0}-\hat{\bm{\Delta}}$ is full-rank.

###### Lemma 6

If $\sigma_{c^{\prime}}({\bm{B}}^{0})>\|\hat{\bm{\Delta}}\|_{2}$ then matrix
${\bm{A}}={\bm{B}}^{0}-\hat{\bm{\Delta}}$ is full-rank.

From Lemma 6 we can see that the success of DPCP-PSGM hinges on how well-
conditioned the matrix ${\bm{B}}^{0}$ is. Specifically, it says that if a
lower-bound on the smallest singular is satisfied then DPCP-PSGM is guaranteed
to converge to the correct complement of the inlier without knowledge of the
correct codimension $c$. From this, we can prove the following Theorem.

###### Theorem 7

Let $\mathbf{B}^{0}\in\mathbb{R}^{D\times c^{\prime}}$ with columns randomly
sampled from a unit $\ell_{2}$ norm spherical distribution where
$c^{\prime}\geq c$ with $c$ denoting the true codimension of the inliers
subspace $\mathcal{S}$ that satisfies Assumption 1. If

$1-C_{1}\sqrt{\frac{c^{\prime}}{D}}-\frac{\epsilon}{\sqrt{D}}>\sqrt{c^{\prime}}\kappa(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}-c_{d})$
(18)

where $\kappa=\max_{i}\frac{M\mu^{0}_{i}}{\beta^{K_{0}/K_{\ast}}(1-r_{i})}$
and
$r_{i}=\frac{\left(1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\right)}{1-\mu^{0}_{i}Mc_{D}}\beta^{1/K_{\ast}}$
then with probability at least $1-2\exp(-\epsilon^{2}C_{2})$ (where
$C_{1},C_{2}$ are absolute constants), Algorithm 1 with a piecewise
geometrically diminishing step size rule will converge to a matrix
$\hat{{\bm{B}}}^{\ast}$ such that
$\operatorname*{span}(\hat{{\bm{B}}}^{\ast})\equiv\mathcal{S}_{\perp}$.

Note that quantities $\beta,K_{\ast},K_{0}$ are used in the step-size update
rule that is used as defined in (8) (See also full version of Theorem 1 in
Appendix)). Theorem 7 shows that we can randomly initialize DPCP-PSGM, with a
matrix $\hat{{\bm{B}}}^{0}$ whose number of columns $c^{\prime}$ is an
overestimate of the true codimension $c$ of the inliers’ subspace and with
columns sampled independently by a uniform distribution over the unit sphere
and recover a matrix that will span the orthogonal complement of
$\mathcal{S}$. The probability of success depends on the geometry of the
problem since condition (18) is trivially satisfied (RHS of (18) tends to 0
since $\bm{\eta}_{\bm{\mathcal{O}}}\rightarrow 0$ and
$\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}\rightarrow c_{d}$) in the
continuous case which incurs a benign geometry. Moreover, the a less benign
geometry would increase the value of
$(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}-c_{d})$
thus requiring a smaller initial codimension $c^{\prime}$ that would lead to
larger values the LHS of (18).

##  5 Numerical Simulations

In this section we demonstrate the effectiveness of the proposed DPCP
formulation and the derived DPCP-PSGM algorithm in recovering orthogonal
complements of subspace of unknown codimension. We compare the proposed
algorithm with previously developed methods i.e., DPCP-IRLS Tsakiris & Vidal
(2018) and the Riemannian Subgradient Method (RSGM) Zhu et al. (2019). Recall
that both DPCP-IRLS and RSGM address DPCP problem by enforcing orthogonality
constraints, and thus both algorithms are quite sensitive if the true
codimension of $\mathcal{S}$ is not known a priori. Further, they are both
initialized using of spectral initialization i.e.,
$\hat{{\bm{B}}}^{0}\in\mathbb{R}^{D\times c^{\prime}}$ which contains the
first $c^{\prime}$ eigenvectors of matrix
$\bm{\tilde{\mathcal{X}}}\bm{\tilde{\mathcal{X}}}^{\top}$ as its columns, as
proposed in Tsakiris & Vidal (2018); Zhu et al. (2019).

Robustness to outliers in the unknown codimension regime. In this experiment
we set the dimension of the ambient space to $D=200$. We randomly generate $N$
inliers uniformly distributed with unit $\ell_{2}$ norm in a $d=195$
dimensional subspace (hence for its codimension we have $c=D-d=5$). Following
a similar process we generate $M$ outliers that live in the ambient space and
are sampled from a uniform distribution over the unit sphere. Fig. 2,
illustrates the distances (see Appendix) of the recovered matrix
$\hat{\bm{B}}$ as obtained by the proposed DPCP-PSGM algorithm initialized
with an overestimate $c^{\prime}=10$ of the true codimension $c$ codimension
and two versions of RSGM i.e., RSGM when it is given as input true $c=5$ and
RSGM when being incognizant of $c$ and hence it initialized with a
$c^{\prime}=10$ of $c$As is shown in Fig. 2 (right), RSGM fails to recover the
correct orthogonal complement of $\mathcal{S}$ when it is provided with an
overestimate of the true $c$ which is attributed to spectral initialization
and the imposed orthogonality constraints. On the contrary, DPCP-PSGM displays
a remarkably robust behavior (Fig. 2(middle)) even without knowing the true
value of $c$, performing similarly to RSGM when the latter knows beforehand
the correct codimension (Fig. 2(left)).

|  |
---|---|---
Figure 2: Distances of the recovered $\hat{\bm{B}}$ from the true orthogonal
complements $\mathcal{S}_{\perp}$ as recovered by the proposed DPCP-PSGM
algorithm provided an overestimated of the true $c$ i.e., $c^{\prime}=10$
(left), RSGM provided the true $c$ (middle) and RSGM provided $c^{\prime}=10$
(right). Darker colors reflect higher values of distances while lighter colors
indicate successful recoveries of ${\bm{B}}$.

Recovery of the true codimension. Here we test DPCP-PSGM on the recovery of
the true codimension $c$ of the inliers’ subspace $\mathcal{S}$. Again, we set
$D=200$ and follow the same process described above for generating $N=1500$
inliers. We vary the true codimension of $\mathcal{S}$ from $c=10$ to $20$ and
consider two different outlier’s ratios $r$, defined as $r=\frac{M}{M+N}$,
namely $r=0.6$ and $r=0.7$. In both cases, DPCP-PSGM is initialized with the
same overestimate of $c$ i.e., $c^{\prime}=30$. In Fig. 3 we report the
estimated codimensions obtained by DPCP-PSGM for 10 independent trials of the
experiments. It can be observed that DPCP-PSGM achieves 100% for all different
codimensions for $r=0.6$. Moreover, it shows a remarkable performance in
estimating the correct $c$’s even in the more challenging case corresponding
to outliers’ ratios equal to 0.7. with the estimated codimensions being close
to the true values even in the cases that it fails to exactly compute $c$. The
results corroborate the theory showing that the DPCP-PSGM with random
initialization biases the solutions of $\hat{\mathbf{B}}$ towards matrices
with rank $c$.

|
---|---
Figure 3: Estimated by DPCP-PSGM codimensions for two different outliers’
ratios $r=\frac{M}{M+N}$ (a) $r=0.6$ (left) and (b) $r=0.7$ (right)

A hyperspectral imaging experiment See Appendix B.1.

##  6 Conclusions

We proposed a simple framework which allows us to perform robust subspace
recovery without requiring a priori knowledge of the subspace codimension.
This is based on Dual Principal Component Pursuit (DPCP) and thus is amenable
to handling subspaces of high relative dimensions. We observed that a
projected subgradient methods (PSGM) induces implicit bias and converges to a
matrix that spans a basis of the orthogonal complement of the inliers subspace
even as long as a) we overestimate it codimension, b) lift orthogonality
constraints enforced in previous DPCP formulations and c) use random
initialization. We provide empirical results that corroborate the developed
theory and showcase the merits of our approach.

#### Ethics Statement

This work focuses on theoretical aspects of robust subspace recovery problem
which is a well-established topic in machine learning research. The research
conducted in the framework of this work raises no ethical issues or any
violations vis-a-vis the ICLR Code of Ethics.

## References

  * Bai et al. (2019) Yu Bai, Qijia Jiang, and Ju Sun.  Subgradient descent learns orthogonal dictionaries.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=HklSf3CqKm.
  * Candès et al. (2011) Emmanuel J Candès, Xiaodong Li, Yi Ma, and John Wright.  Robust principal component analysis?  _Journal of the ACM (JACM)_ , 58(3):1–37, 2011\.
  * Ding et al. (2021) Tianyu Ding, Zhihui Zhu, Rene Vidal, and Daniel P Robinson.  Dual principal component pursuit for robust subspace learning: Theory and algorithms for a holistic approach.  In Marina Meila and Tong Zhang (eds.), _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pp. 2739–2748. PMLR, 18–24 Jul 2021.  URL https://proceedings.mlr.press/v139/ding21b.html.
  * Fischler & Bolles (1981) Martin A Fischler and Robert C Bolles.  Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography.  _Communications of the ACM_ , 24(6):381–395, 1981\.
  * Giampouras et al. (2019) Paris V. Giampouras, Athanasios A. Rontogiannis, and Konstantinos D. Koutroumbas.  Alternating iteratively reweighted least squares minimization for low-rank matrix factorization.  _IEEE Transactions on Signal Processing_ , 67(2):490–503, 2019.  doi: 10.1109/TSP.2018.2883921.
  * Gunasekar et al. (2018) Suriya Gunasekar, Jason Lee, Daniel Soudry, and Nathan Srebro.  Characterizing implicit bias in terms of optimization geometry.  In Jennifer Dy and Andreas Krause (eds.), _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 1832–1841. PMLR, 10–15 Jul 2018.  URL https://proceedings.mlr.press/v80/gunasekar18a.html.
  * Jolliffe & Cadima (2016) Ian T Jolliffe and Jorge Cadima.  Principal component analysis: a review and recent developments.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 374(2065):20150202, 2016.
  * Lerman & Maunu (2018) Gilad Lerman and Tyler Maunu.  An overview of robust subspace recovery.  _Proceedings of the IEEE_ , 106(8):1380–1410, 2018.
  * Lerman et al. (2015) Gilad Lerman, Michael B McCoy, Joel A Tropp, and Teng Zhang.  Robust computation of linear models by convex relaxation.  _Foundations of Computational Mathematics_ , 15(2):363–410, 2015.
  * Maunu et al. (2019) Tyler Maunu, Teng Zhang, and Gilad Lerman.  A well-tempered landscape for non-convex robust subspace recovery.  _Journal of Machine Learning Research_ , 20(37), 2019.
  * McCoy & Tropp (2011) Michael McCoy and Joel A Tropp.  Two proposals for robust PCA using semidefinite programming.  _Electronic Journal of Statistics_ , 5:1123–1160, 2011\.
  * Qu et al. (2020) Qing Qu, Zhihui Zhu, Xiao Li, Manolis C Tsakiris, John Wright, and René Vidal.  Finding the sparsest vectors in a subspace: Theory, algorithms, and applications.  _arXiv preprint arXiv:2001.06970_ , 2020.
  * Rahmani & Atia (2017) Mostafa Rahmani and George Atia.  Coherence pursuit: Fast, simple, and robust subspace recovery.  In Doina Precup and Yee Whye Teh (eds.), _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pp. 2864–2873. PMLR, 06–11 Aug 2017.  URL https://proceedings.mlr.press/v70/rahmani17a.html.
  * Tsakiris & Vidal (2018) Manolis C. Tsakiris and René Vidal.  Dual principal component pursuit.  _Journal of Machine Learning Research_ , 19(18):1–50, 2018.
  * Vershynin (2010) Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  _arXiv preprint arXiv:1011.3027_ , 2010.
  * Vidal et al. (2016) René Vidal, Yi Ma, and S Shankar Sastry.  _Generalized principal component analysis_ , volume 5.  Springer, 2016.
  * Xu et al. (2012) Huan Xu, Constantine Caramanis, and Sujay Sanghavi.  Robust PCA via outlier pursuit.  _IEEE Transactions on Information Theory_ , 58(5):3047–3064, 2012.
  * You et al. (2017a) C. You, D. Robinson, and R. Vidal.  Provable self-representation based outlier detection in a union of subspaces.  In _IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 4323–4332, 2017a.
  * You et al. (2017b) Chong You, Daniel P Robinson, and René Vidal.  Provable self-representation based outlier detection in a union of subspaces.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 3395–3404, 2017b.
  * You et al. (2020) Chong You, Zhihui Zhu, Qing Qu, and Yi Ma.  Robust recovery via implicit bias of discrepant learning rates for double over-parameterization.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin (eds.), _Advances in Neural Information Processing Systems_ , volume 33, pp. 17733–17744. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/cd42c963390a9cd025d007dacfa99351-Paper.pdf.
  * Zhang & Lerman (2014) Teng Zhang and Gilad Lerman.  A novel m-estimator for robust PCA.  _The Journal of Machine Learning Research_ , 15(1):749–808, 2014.
  * Zhu et al. (2019) Z. Zhu, T. Ding, M. C. Tsakiris, D. P. Robinson, and R. Vidal.  A linearly convergent method for non-smooth non-convex optimization on the grassmannian with applications to robust subspace and dictionary learning.  In _Neural Information Processing Systems (NIPS)_ , 2019.
  * Zhu et al. (2018) Zhihui Zhu, Yifan Wang, Daniel Robinson, Daniel Naiman, René Vidal, and Manolis Tsakiris.  Dual principal component pursuit: Improved analysis and efficient algorithms.  In _Advances in Neural Information Processing Systems 2018_ , volume 31. Curran Associates, Inc., 2018.  URL https://proceedings.neurips.cc/paper/2018/file/af21d0c97db2e27e13572cbf59eb343d-Paper.pdf.

##  Appendix A Appendix

###### Theorem 1

(Theorem 3 of Zhu et al. (2018)) Let $\\{\hat{{\bm{b}}}_{k}\\}$ the sequence
generated by the projected subgradient method in Zhu et al. (2018), with
initialization $\hat{{\bm{b}}}_{0}$ such that

$\theta_{0}<\mathrm{arctan}\left(\frac{Nc_{{\bm{X}},\min}}{N\eta_{{\bm{X}}}+M\eta_{{\bm{O}}}}\right)$
(19)

where $\theta_{0}$ denotes the principal angle of ${\bm{b}}^{0}$ from
$\mathcal{S}_{\perp}$, and

$Nc_{{\bm{X}},\min}\geq N\eta_{{\bm{X}}}+M\eta_{{\bm{O}}}$ (20)

Let
$\mu^{\prime}\coloneqq\frac{1}{4\max\\{Nc_{{\bm{X}},\min},Mc_{{\bm{O}},\max}\\}}$.
If $\mu^{0}\leq\mu^{\prime}$ and the step size $\mu^{k}$ is updated according
to a piece-wise geometrically diminishing rule given as

$\mu^{k}=\left\\{\begin{array}[]{l l}\mu^{0},&k<K_{0}\\\
\mu^{0}\beta^{\lfloor(k-K_{0})/K_{\ast}\rfloor+1},&k\geq
K_{0}\end{array}\right.$ (21)

where $\beta<1$, $\lfloor\cdot\rfloor$ is the floor function, and
$K_{0},K_{\ast}\in\mathbb{N}$ are chosen such that

$\begin{split}&K_{0}\geq K^{\diamondsuit}(\mu^{0}),\\\
&K_{\ast}\geq\left(\sqrt{2}\beta\mu^{\prime}\left(Nc_{{\bm{X}}}-\left(N\eta_{{\bm{X}}}+M\eta_{{\bm{O}}}\right)\right)\right)^{-1}\end{split}$

where,

$K^{\diamondsuit}(\mu)\coloneqq\frac{\tan(\theta_{0})}{\mu\left(Nc_{{\bm{X}},\min}-\max\\{1,\tan(\theta_{0}\\}\left(N\eta_{{\bm{X}}}+M\eta_{{\bm{O}}}\right)\right)}$
(22)

then for the angle $\theta_{k}$ between $\hat{{\bm{b}}}^{k}$ and
$\mathcal{S}_{\perp}$ it holds

$\tan(\theta_{0})\leq\left\\{\begin{array}[]{l
l}\max\\{\tan(\theta_{0}),\frac{\mu_{0}}{\sqrt{2}\mu^{\prime}}\\},&k<K_{0}\\\
\frac{\mu^{0}}{\sqrt{2}\mu^{\prime}}\beta^{\lfloor(k-K_{0})/K_{\ast}\rfloor},&k\geq
K_{0}\end{array}\right.$ (23)

###  A.1 Proof of Lemma 2

###### Lemma 2

In the continuous case, the discrete DPCP problem given in (3) is reformulated
as,

$\begin{split}\min_{{\bm{B}}\in\mathbb{R}^{D\times
c^{\prime}}}\sum^{c^{\prime}}_{i=1}\left(p\mathbb{E}_{\bm{\mu}_{\mathbb{S}^{D-1}}}[f_{\mathbf{b}_{i}}]+(1-p)\mathbb{E}_{\bm{\mu}_{\mathbb{S}^{D-1}\cap\bm{S}}}[f_{{\bm{b}}_{i}}]\right)=\sum^{c^{\prime}}_{i=1}\|{\bm{b}}_{i}\|_{2}\left(pc_{D}+(1-p)c_{d}\mathrm{cos}(\phi_{i})\right)\\\
\ \ \mathrm{s.t.}\ \
\|{\bm{b}}_{i}\|_{2}=1,~{}~{}i=1,2,\dots,c^{\prime}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(24)

where $f_{{\bm{b}}}:\mathbb{S}^{D-1}\rightarrow\mathbb{R}$,
$f_{{\bm{b}}}(\mathbf{z})=|\mathbf{z}^{\top}{\bm{b}}|$, $\phi_{i}$ is the
principal angle of ${\bm{b}}_{i}$ from the inliers subspace $\mathcal{S}$ and
$p$ is the probability of occurrence of an outlier.

Proof  We define the discrete measures
$\bm{\mu}_{\bm{\mathcal{X}}},\bm{\mu}_{\bm{\mathcal{O}}}$ associated with the
inliers and outliers, respectively as,

$\bm{\mu}_{\bm{\mathcal{X}}}(\mathbf{z})=\frac{1}{N}\sum^{N}_{j=1}\delta(\mathbf{z}-\bm{\mathit{o}}_{j}),~{}~{}\bm{\mu}_{\bm{\mathcal{O}}}(\mathbf{z})=\frac{1}{M}\sum^{N}_{j}\delta(\mathbf{z}-\bm{\mathit{x}}_{j})$
(25)

where $\delta(\cdot)$ is the Dirac function. Recall that,

$\int_{\mathbf{z}\in\mathbb{S}^{D-1}}g(\mathbf{z})\delta(\mathbf{z}-\mathbf{z}_{0})d\bm{\mu}_{\mathbb{S}^{D-1}}=g(z_{0})$
(26)

where $g:\mathbb{S}^{D-1}\rightarrow\mathbb{R}$ and
$\bm{\mu}_{\mathbb{S}^{D-1}}$ is the uniform measure on $\mathbb{S}^{D-1}$.

The DPCP objective for the discrete version of the problem divided by $M+N$
can be written as,

$\begin{split}&\frac{1}{M+N}\sum^{c^{\prime}}_{i=1}\|\bm{\tilde{\mathcal{X}}}^{\top}\mathbf{b}_{i}\|_{1}=\frac{1}{M+N}\sum^{c^{\prime}}_{i=1}\left(\|\bm{\mathcal{X}}^{\top}\mathbf{b}_{i}\|_{1}+\|\bm{\mathcal{O}}^{\top}\mathbf{b}_{i}\|_{1}\right)=\frac{1}{M+N}\sum^{c^{\prime}}_{i}\left(\sum^{N}_{j=1}\
|\bm{\mathit{x}}^{\top}_{j}\mathbf{b}_{i}|+\sum^{M}_{j=1}|\bm{\mathit{o}}^{\top}_{j}\mathbf{b}_{i}|\right)\\\
&=\frac{1}{M+N}\sum^{c^{\prime}}_{i=1}\left(\sum^{N}_{j=1}\int_{\mathbf{z}\in\mathbb{S}^{D-1}}|\mathbf{z}^{\top}\mathbf{b}_{i}|\delta(\mathbf{z}-\bm{\mathit{x}}_{j})d\bm{\mu}_{\mathbb{S}^{D-1}}+\sum^{M}_{j}\int_{\mathbf{z}\in\mathbb{S}^{D-1}}|\mathbf{z}^{\top}\mathbf{b}_{i}|\delta(\mathbf{z}-\bm{\mathit{o}}_{j})d\bm{\mu}_{\mathbb{S}^{D-1}}\right)\\\
&=\frac{1}{M+N}\sum^{c^{\prime}}_{i=1}\left(\int_{\mathbf{z}\in\mathbb{S}^{D-1}}|\mathbf{z}^{\top}\mathbf{b}_{i}|\sum^{N}_{j=1}\delta(\mathbf{z}-\bm{\mathit{x}}_{j})d\bm{\mu}_{\mathbb{S}^{D-1}}+\int_{\mathbf{z}\in\mathbb{S}^{D-1}}|\mathbf{z}^{\top}\mathbf{b}_{i}|\sum^{M}_{j}\delta(\mathbf{z}-\bm{\mathit{o}}_{j})d\bm{\mu}_{\mathbb{S}^{D-1}}\right)\\\
&=\sum^{c^{\prime}}_{i=1}\left(p\mathbb{E}_{\bm{\mu}_{\bm{\mathcal{X}}}}[f_{\mathbf{b}_{i}}]+(1-p)\mathbb{E}_{\bm{\mu}_{\bm{\mathcal{O}}}}[f_{\mathbf{b}_{i}}]\right)\end{split}$

Note that $\bm{\mu}_{\bm{\mathcal{X}}},\bm{\mu}_{\bm{\mathcal{O}}}$ arise by
discretizing the continuous uniform measures $\bm{\mu}_{\mathbb{S}^{D-1}}$ and
$\bm{\mu}_{\mathbb{S}^{D-1}\cap\bm{S}}$ respectively
($\bm{\mu}_{\mathbb{S}^{D-1}\cap\bm{S}}$ denotes the uniform measure on
$\mathbb{S}^{D-1}\cap\mathcal{S}$) and $p$ is the probability of occurrence of
an outlier i.e., $\frac{M}{M+N}\rightarrow p$ as $M,N\rightarrow\infty$ (
$1-p$ corresponds to the probability of occurrence of an inlier). That being
said, the continuous version of DPCP can be simply stated by replacing
$\bm{\mu}_{\bm{\mathcal{X}}},\bm{\mu}_{\bm{\mathcal{O}}}$ with
$\bm{\mu}_{\mathbb{S}^{D-1}}$ and $\bm{\mu}_{\mathbb{S}^{D-1}\cap\bm{S}}$ in
equation A.1 as follows,

$\min_{\mathbf{B}}\sum^{c^{\prime}}_{i=1}\left(p\mathbb{E}_{\bm{\mu}_{\mathbb{S}^{D-1}}}[f_{\mathbf{b}_{i}}]+(1-p)\mathbb{E}_{\bm{\mu}_{\mathbb{S}^{D-1}\cap\bm{S}}}[f_{\mathbf{b}_{i}}]\right)$
(27)

The RHS of (9) immediately shows up by invoking Proposition 4 in Tsakiris &
Vidal (2018).

###  A.2 Proof of Lemma 3

Lemma 3: A projected subgradient algorithm consisting of the steps described
in (10) using a piecewise geometrically diminishing step size rule (see (8) in
Theorem 1) will almost surely asymptotically converge to a matrix
$\hat{{{\bm{B}}}}^{\ast}\in\mathbb{R}^{D\times c^{\prime}}$ whose columns
$\hat{{\bm{b}}}^{\ast}_{i},~{}i=1,2,\dots,c^{\prime}$ will be normal vectors
of the inliers’ subspace when randomly initialized with vectors
${\bm{b}}^{0}_{i}\in\mathbb{S}^{D-1},~{}i=1,2,\dots,c^{\prime}$ uniformly
distributed over the sphere $\mathbb{S}^{D-1}$.
Proof:

The proof can be trivially obtained by noticing a) that the condition for
convergence i.e., inequality (7) of the projected subgradient algorithm given
in Theorem 1 becomes $\theta^{0}_{i}<\frac{\pi}{2}$ in the continuous case
(since $\eta_{{\bm{X}}}\rightarrow 0$, $\eta_{{\bm{O}}}\rightarrow 0$,
$c_{{\bm{X}},\min}\rightarrow c_{d}>0$) and b) the set of unit $\ell_{2}$-norm
vectors ${\bm{b}}^{0}_{i}$s, $i=1,2,\dots,c^{\prime}$ sampled independently by
a uniform distribution over the sphere and whose principal angle
$\theta^{0}_{i}$ is $\frac{\pi}{2}$ form the inliers’ subspace has measure
$0$. $\blacksquare$

###  A.3 Proof of Theorem 5

###### Lemma 4

The PSGM iterates
$\hat{{\bm{b}}}^{k}_{i},~{}i=1,2,\dots,c^{\prime},~{}k=1,2,\dots$ given in
(10), when randomly initialized with $\hat{{\bm{b}}}^{0}_{i}$s,
$i=1,2,\dots,c^{\prime}$ that are independently drawn from a spherical
distribution with unit $\ell_{2}$ norm converge almost surely to $c^{\prime}$
normal vectors of the inliers subspace $\mathcal{S}$ denoted as
$\hat{{\bm{b}}}^{\ast}_{i},~{}i=1,2,\dots,c^{\prime}$ that are given by

$\hat{{\bm{b}}}^{\ast}_{i}=\frac{\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{i})}{\|\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{i})\|_{2}},~{}~{}~{}~{}~{}i=1,2,\dots,c^{\prime}$
(28)

Proof  Let us assume ${\bm{b}}^{0}_{i}=\hat{{\bm{b}}}^{0}_{i}$. The iterates
of subgradients steps of PSGM can be written in the following form,

$\begin{split}{\bm{b}}^{1}_{i}&=(1-\mu^{0}_{i}pc_{D})\hat{{\bm{b}}}^{0}_{i}-\mu^{0}_{i}(1-p)c_{d}\hat{\mathbf{s}}\\\
{\bm{b}}^{2}_{i}&=(1-\mu^{1}_{i}pc_{D})\hat{{\bm{b}}}^{1}_{i}-\mu^{1}_{i}(1-p)c_{d}\hat{\mathbf{s}}\\\
\vdots~{}~{}&=~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\vdots\\\
{\bm{b}}^{K}_{i}&=(1-\mu^{K-1}_{i}pc_{D})\hat{{\bm{b}}}^{K-1}_{i}-\mu^{K-1}_{i}(1-p)c_{d}\hat{\mathbf{s}}\end{split}$
(29)

By projecting each update of PSGM onto $\mathcal{S}_{\perp}$ and since
$\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{k}_{i})=\|{\bm{b}}^{k}_{i}\|\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{k}_{i})$
we have,

$\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{k+1}_{i})=\frac{(1-\mu^{k}_{i}pc_{D})}{\|{\bm{b}}^{k}_{i}\|}\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{k}_{i})$
(30)

We can thus easily derive the following form for
$\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{K}_{i})$,

$\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{K}_{i})=\left(\prod^{K-1}_{k=1}\frac{(1-\mu^{k}_{i}pc_{D})}{\|{\bm{b}}^{k}_{i}\|_{2}}\right)\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{b}}^{0}_{i})$
(31)

We know from Theorem 1 and Lemma 3 when DPCP-PSGM is initialized with
${\bm{b}}^{0}_{i},~{}i=1,2,\dots,c^{\prime}$s randomly drawn according
according to a spherical distribution then it will almost surely converge as
$K\rightarrow\infty$ to vectors
$\hat{{\bm{b}}}^{\ast}_{i},~{}i=1,2,\dots,c^{\prime}$ i.e.,
$\hat{{\bm{b}}}^{K}_{i}\rightarrow\hat{{\bm{b}}}^{\ast}_{i}$ where
$\hat{{\bm{b}}}^{\ast}_{i}\in\mathcal{S}_{\perp}$. Hence
$\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{K}_{i})\rightarrow\hat{{\bm{b}}}^{\ast}_{i}$
as $K\rightarrow\infty$. Note that from Theorem 1 we have that
$\mu^{k}_{i}\neq\frac{1}{pc_{D}}~{}~{}~{}\forall k=\\{1,2,\dots,K\\}$ hence
$\prod^{K-1}_{k=1}\frac{(1-\mu^{k}_{i}pc_{D})}{\|\hat{{\bm{b}}}^{k}_{i}\|_{2}}\neq
0$. From 31 and after projecting on the unit sphere and we thus have
$\hat{{\bm{b}}}^{\ast}_{i}=\frac{\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{i})}{\|\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{i})\|_{2}}$.

###### Theorem 5

Let $\hat{{\bm{B}}^{0}}\in\mathbb{R}^{D\times c^{\prime}}$ where
$c^{\prime}\geq c$ with $c$ denoting the true codimension of the inliers
subspace $\mathcal{S}$, consisting of unit $\ell_{2}$ norm column vectors
$\hat{{\bm{b}}}^{0}_{i}\in\mathbb{S}^{D-1},i=1,2,\dots,c^{\prime}$ that are
independently drawn from uniform distribution over the sphere
$\mathbb{S}^{D-1}$. A PSGM algorithm initialized with $\hat{{\bm{B}}^{0}}$
will almost surely converge to a matrix $\hat{{\bm{B}}}^{\ast}$ such that
$\operatorname*{span}(\hat{{\bm{B}}}^{\ast})\equiv\mathcal{S}_{\perp}$.

Proof  From Lemma 4 we have that for each initial unit norm vector
${\bm{b}}^{0}_{i}$ which corresponds to the $i$th column of ${\bm{B}}^{0}$
will almost surely converge to
$\hat{{\bm{b}}}^{\ast}_{i}=\frac{\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{i})}{\|\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{i})\|_{2}}$.
We can thus write
${\bm{B}}^{\ast}=\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{B}}^{0})\bm{\Gamma}$
where $\bm{\Gamma}$ is a full-rank diagonal matrix given be
$\bm{\Gamma}=\operatorname*{diag}(\frac{1}{\|\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{1})\|_{2}},\frac{2}{\|\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{2})\|_{2}},\dots,\frac{1}{\|\mathcal{P}_{\mathcal{S}_{\perp}}(\hat{{\bm{b}}}^{0}_{c^{\prime}})\|_{2}})$.
Note that $\mathcal{P}_{\mathcal{S}_{\perp}}$ is a linear projection and thus
we can write
$\mathcal{P}_{\mathcal{S}_{\perp}}({\bm{B}}^{0})={\bm{B}}_{\mathcal{S}_{\perp}}{\bm{B}}^{\top}_{\mathcal{S}_{\perp}}$
where $\mathbf{B}_{\mathcal{S}_{\perp}}\mathbb{R}^{D\times c}$ is an
orthonormal matrix which spans $\mathcal{S}_{\perp}$. Note that the
probability of sampling a low-rank matrix
${\bm{B}}_{0}=[{\bm{b}}^{0}_{1},{\bm{b}}^{0}_{2},\dots,{\bm{b}}^{0}_{c^{\prime}}]$
when columns ${\bm{b}}^{0}_{i}$s are randomly and independently drawn from a
spherical distribution is zero. We thus have
${\bm{B}}^{\ast}={\bm{B}}_{\mathcal{S}_{\perp}}{\bm{B}}^{\top}_{\mathcal{S}_{\perp}}{\bm{B}}^{0}\bm{\Gamma}$
with $\operatorname*{rank}({\bm{B}}^{\ast})=c$.

###### Lemma 6

If $\sigma_{c^{\prime}}({\bm{B}}^{0})>\|\hat{\bm{\Delta}}\|_{2}$ then the rank
of matrix ${\bm{B}}^{0}-\hat{\bm{\Delta}}$ equals $c^{\prime}$.

Proof  Let ${\bm{A}}={\bm{B}}^{0}-\hat{\bm{\Delta}}$. From singular value
perturbation inequalities we have
$|\sigma_{i}({\bm{A}})-\sigma_{i}({\bm{B}}^{0})|\leq\|\hat{\bm{\Delta}}\|_{2}$,
for $i=1,2,\dots,c^{\prime}$. Hence it holds,

$-\|\hat{\bm{\Delta}}\|_{2}\leq\sigma_{i}({\bm{A}})-\sigma_{i}({\bm{B}}^{0})$
(32)

If $\sigma_{c^{\prime}}({\bm{B}}^{0})>\|\hat{\bm{\Delta}}\|_{2}$ then from
equation 32 we get

$\begin{split}\sigma_{c^{\prime}}({\bm{A}})>0\end{split}$ (33)

hence the matrix ${\bm{B}}^{0}-\hat{\bm{\Delta}}$ will be full-rank.

###  A.4 Proof of Theorem 7

We first give the following Lemmas:

###### Lemma 7

For the $\ell_{2}$ norm of $\bm{\mathit{e}}_{{\bm{O}}}^{i,k}$ for any
$k=1,2,\dots,K$ and $i=1,2,\dots,c^{\prime}$ it holds,

$\|\bm{\mathit{e}}_{{\bm{O}}}^{i,k}\|_{2}\leq\eta_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}-c_{d}$
(34)

Proof

$\begin{split}\|\bm{\mathit{e}}_{{\bm{O}}}^{i,k}\|_{2}=\|\bm{\mathit{o}}_{b}-c_{d}\mathbf{b}\|_{2}&=\|\frac{1}{M}\bm{\mathcal{O}}\mathrm{Sgn}(\bm{\mathcal{O}}^{\top}{\bm{b}})-c_{d}{\bm{b}}\|_{2}\\\
&=\|\frac{1}{M}(\mathbf{I}-{\bm{b}}{\bm{b}}^{\top})\bm{\mathcal{O}}\mathrm{Sgn}(\bm{\mathcal{O}}^{\top}{\bm{b}})+\frac{1}{M}\mathbf{b}\mathbf{b}^{\top}\bm{\mathcal{O}}\mathrm{Sgn}(\bm{\mathcal{O}}^{\top}{\bm{b}})-c_{d}{\bm{b}}\|_{2}\\\
&\leq\|\frac{1}{M}(\mathbf{I}-{\bm{b}}{\bm{b}}^{\top})\bm{\mathcal{O}}\mathrm{Sgn}(\bm{\mathcal{O}}^{\top}{\bm{b}})\|_{2}+(\frac{1}{M}\|\bm{\mathcal{O}}^{\top}{\bm{b}}\|_{1}-c_{d})\|{\bm{b}}\|_{2}\\\
&\leq\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}-c_{d}\end{split}$
(35)

where we have used the fact that $\|{\bm{b}}\|_{2}=1$.

###### Lemma 8

Let the step size of Algorithm 1 (DPCP-PSGM) $\mu^{k}_{i}$ being updated
following the piecewise geometrically diminishing step size rule with

$\beta<\left(\frac{1-\mu^{0}_{i}Mc_{D}}{1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)}\right)^{K_{\ast}}.$

For the spectral norm of $\hat{\bm{\Delta}}$ it holds
$\|\hat{\bm{\Delta}}\|_{2}<\sqrt{c^{\prime}}\kappa(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}-c_{d})$
where $\kappa=\max_{i}\frac{M\mu^{0}_{i}}{\beta^{K_{0}/K_{\ast}}(1-r_{i})}$
and
$r_{i}=\frac{\left(1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\right)}{1-\mu^{0}_{i}Mc_{D}}\beta^{1/K_{\ast}}$.

Proof  We first bound the $\ell_{2}$ norm of vectors $\mathbf{b}^{j}_{i}$s. We
have that $\forall i=1,2,\dots,c^{\prime}$ and $j=1,2,\dots,K$ it holds

${\bm{b}}^{j+1}_{i}=\hat{{\bm{b}}}^{j}_{i}-\mu^{j}_{i}\left(\bm{\mathcal{X}}\mathrm{Sgn}(\bm{\mathcal{X}}^{\top}\hat{{\bm{b}}}^{j}_{i})+\bm{\mathcal{O}}\mathrm{Sgn}(\bm{\mathcal{O}}^{\top}\hat{{\bm{b}}}^{j}_{i})\right)$
(36)

We define the quantities

$\bm{\eta}_{\bm{\mathcal{X}}}\coloneqq\max_{\mathbf{b}\in\mathbb{S}^{D-1}}\frac{1}{N}\|(\mathcal{P}_{\mathcal{S}}-\hat{\mathbf{b}}^{j}_{i}\hat{\mathbf{b}}^{j,\top}_{i})\bm{\mathcal{X}}\mathrm{Sgn}(\bm{\mathcal{X}}^{\top}\hat{\mathbf{b}}^{j}_{i})\|_{2}$
(37)
$c_{{\bm{X}},\max}\coloneqq\max_{\mathbf{b}\in\mathbb{S}^{D-1}}\frac{1}{N}\|\bm{\mathcal{X}}^{\top}\hat{\mathbf{b}}^{j}_{i}\|_{1}$
(38)

$\|\hat{\mathbf{b}}^{j}_{i}\|_{2}=1$ hence

$\begin{split}\|{\bm{b}}^{j+1}_{i}\|_{2}&\leq
1+\mu^{j}_{i}\|\bm{\mathcal{X}}\mathrm{Sgn}(\bm{\mathcal{X}}^{\top}\hat{{\bm{b}}}^{j}_{i})+\bm{\mathcal{O}}\mathrm{Sgn}(\bm{\mathcal{O}}^{\top}\hat{{\bm{b}}}^{j}_{i}\|_{2}\\\
&\leq
1+\mu^{j}_{i}\left(\|\bm{\mathcal{X}}\mathrm{Sgn}(\bm{\mathcal{X}}^{\top}\hat{{\bm{b}}}^{j}_{i})\|_{2}+\|\bm{\mathcal{O}}\mathrm{Sgn}(\bm{\mathcal{O}}^{\top}\hat{\mathbf{b}}^{j}_{i})\|_{2}\right)\\\
&\leq
1+\mu^{j}_{i}(\|(\mathcal{P}_{\mathcal{S}}-\hat{\mathbf{b}}^{j}_{i}\hat{{\bm{b}}}^{j,\top}_{i})\bm{\mathcal{X}}\mathrm{Sgn}(\bm{\mathcal{X}}^{\top}\hat{{\bm{b}}}^{j}_{i})\|_{2}+\|(\hat{\mathbf{b}}^{j}_{i}\hat{{\bm{b}}}^{j,\top}_{i})\bm{\mathcal{X}}\mathrm{Sgn}(\bm{\mathcal{X}}^{\top}\hat{{\bm{b}}}^{j}_{i})\|_{2}\\\
&+\|\left(1-\hat{{\bm{b}}}^{j}_{i}\hat{\mathbf{b}}^{j,\top}_{i}\right)\bm{\mathcal{O}}\mathrm{Sgn}(\bm{\mathcal{O}}^{\top}\hat{\mathbf{b}}^{j}_{i})\|_{2}+\|\hat{{\bm{b}}}^{j}_{i}\hat{\mathbf{b}}^{j,\top}_{i}\bm{\mathcal{O}}\mathrm{Sgn}(\bm{\mathcal{O}}^{\top}\hat{{\bm{b}}}^{j}_{i})\|_{2})\\\
&\leq
1+\mu^{j}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\end{split}$

Due to equation A.4 and since $\mu^{j}_{i}$ follows a non-increasing path as
$j\rightarrow k$, the scalar term
$\prod^{k}_{j=0}\frac{\|\mathbf{b}^{j}_{i}\|_{2}}{(1-\mu^{j}_{i}Mc_{D})}$ is
bounded above as follows,

$\prod^{k}_{j=0}\frac{\|{\bm{b}}^{j}_{i}\|_{2}}{(1-\mu^{j}_{i}Mc_{D})}\leq\left(\frac{\left(1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\right)}{1-\mu^{0}_{i}Mc_{D}}\right)^{k}$
(39)

We now focus on the geometrically diminishing step size rule given in equation
8. We have $\mu^{k}_{i}=\mu^{0}_{i}\beta^{\lfloor
k-K_{0}/K_{\ast}\rfloor+1}<\mu^{0}_{i}\beta^{(k-K_{0})/K_{\ast}}$ for $k\geq
K_{0}$ and $\mu^{0}_{i}\beta^{(k-K_{0})/K_{\ast}}>\mu^{0}_{i}$ for $k<K_{0}$.
Hence we can get the following upper bound

$\begin{split}&\lim_{K\rightarrow\infty}\sum^{K-1}_{k=0}\prod^{k}_{j=0}\frac{\|{\bm{b}}^{j}_{i}\|_{2}}{(1-\mu^{j}_{i}Mc_{D})}\mu^{k}_{i}M<\\\
&M\lim_{K\rightarrow\infty}\sum^{K-1}_{k=0}\left(\frac{\left(1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\right)}{1-\mu^{0}_{i}Mc_{D}}\right)^{k}\mu^{0}_{i}\beta^{(k-K_{0})/K_{\ast}}\\\
&\equiv\lim_{K\rightarrow\infty}\sum^{K-1}_{k=0}M\frac{1}{\beta^{K_{0}/K_{\ast}}}\mu^{0}_{i}\left(\frac{\left(1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\right)}{1-\mu^{0}_{i}Mc_{D}}\beta^{1/K_{\ast}}\right)^{k}\end{split}$

The series
$S=\sum^{K-1}_{k=0}\left(\frac{\left(1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\right)}{1-\mu^{0}_{i}Mc_{D}}\beta^{1/K_{\ast}}\right)^{k}$
is geometric and if

$\beta<\left(\frac{1-\mu^{0}_{i}Mc_{D}}{\left(1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\right)}\right)^{K_{\ast}}$
(40)

it converges as $K\rightarrow\infty$ to $\frac{1}{1-r_{i}}$ where
$r_{i}=\frac{\left(1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\right)}{1-\mu^{0}_{i}Mc_{D}}\beta^{1/K_{\ast}}$.

Let us now bound the $\ell_{2}$ norms of the columns of $\hat{\bm{\Delta}}$.
From Lemma 34 we have
$\|\bm{\mathit{e}}_{{\bm{O}}}^{i,k}\|_{2}\leq\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}-c_{d}$.
We can easily thus derive that
$\|\hat{\bm{\delta}}_{i}\|_{2}\leq\kappa(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}-c_{d})$.
For the spectral norm of $\hat{\bm{\Delta}}$ we thus have

$\begin{split}\|\hat{\bm{\Delta}}\|_{2}&=\underset{\bm{\mathit{x}}\neq\mathbf{0}}{\mathrm{sup}}\frac{\|\hat{\bm{\Delta}}\bm{\mathit{x}}\|_{2}}{\|\bm{\mathit{x}}\|_{2}}=\underset{\bm{\mathit{x}}\neq\mathbf{0}}{\mathrm{sup}}\frac{\|\sum^{c^{\prime}}_{i=1}\hat{\bm{\delta}}_{i}x_{i}\|_{2}}{\|\bm{\mathit{x}}\|_{2}}\leq\underset{\bm{\mathit{x}}\neq\mathbf{0}}{\mathrm{sup}}\frac{\sum^{c^{\prime}}_{i=1}\|\hat{\bm{\delta}}_{i}\|_{2}|x_{i}|}{\|\bm{\mathit{x}}\|_{2}}\\\
&\leq\max_{i}\|\hat{\bm{\delta}}_{i}\|_{2}\underset{\bm{\mathit{x}}\neq\mathbf{0}}{\mathrm{sup}}\frac{\|\bm{\mathit{x}}\|_{1}}{\|\bm{\mathit{x}}\|_{2}}\leq\kappa(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}-c_{d})\sqrt{c^{\prime}}\end{split}$
(41)

Where the last inequality arises since
$\|\bm{\mathit{x}}\|_{1}\leq\sqrt{c^{\prime}}\|\bm{\mathit{x}}\|_{2}$.

We then give the Theorem.

###### Theorem 9

(Theorem 5.58 of Vershynin (2010)) Let ${\bm{B}}$ be a $D\times d$ matrix
($D\geq d$) whose columns ${\bm{b}}_{i}$ are independent sub-gaussian
isotropic random vectors in $\mathbb{R}^{D}$ with
$\|{\bm{b}}_{i}\|_{2}=\sqrt{D}$ almost surely. Then for every $t\geq 0$ the
inequality

$\sqrt{D}-C\sqrt{d}-t\leq\sigma_{\min}({\bm{B}})\leq\sigma_{\max}({\bm{B}})\leq\sqrt{D}+C\sqrt{d}+t$
(42)

with probability at least $1-2\exp(-ct^{2})$, where
$C=C^{\prime}_{k},c=c^{\prime}_{K}>0$ depend only on the subgaussian norm
$K=\max_{j}\|{\bm{b}}_{i}\|_{\psi_{2}}$ of the columns.

The proof of Theorem 7 follows next.

###### Theorem 7

Let $\mathbf{B}^{0}\in\mathbb{R}^{D\times c^{\prime}}$ with columns randomly
sampled from a unit $\ell_{2}$ norm spherical distribution where
$c^{\prime}\geq c$ with $c$ denoting the true codimension of the inliers
subspace $\mathcal{S}$ that satisfies Assumption 1. If

$1-C_{1}\sqrt{\frac{c^{\prime}}{D}}-\frac{\epsilon}{\sqrt{D}}>\sqrt{c^{\prime}}\kappa(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}}-c_{d})$
(43)

where $\kappa=\max_{i}\frac{M\mu^{0}_{i}}{\beta^{K_{0}/K_{\ast}}(1-r_{i})}$
and
$r_{i}=\frac{\left(1+\mu^{0}_{i}\left(N(\bm{\eta}_{\bm{\mathcal{X}}}+c_{\bm{\mathcal{X}},\max})+M(\bm{\eta}_{\bm{\mathcal{O}}}+\mathit{c}_{\bm{\mathcal{O}},\mathrm{max}})\right)\right)}{1-\mu^{0}_{i}Mc_{D}}\beta^{1/K_{\ast}}$
then with probability at least $1-2\exp(-\epsilon^{2}C_{2})$ (where
$C_{1},C_{2}$ are absolute constants), Algorithm 1 with a geometrically
diminishing step size rule will converge to a matrix $\hat{{\bm{B}}}^{\ast}$
such that
$\operatorname*{span}(\hat{{\bm{B}}}^{\ast})\equiv\mathcal{S}_{\perp}$.

Proof  By Assumption 1 we have that all columns of ${\bm{B}}_{0}$ will satisfy
the sufficient condition for convergence of DPCP-PSGM (Algorithm 1) to a
normal vector of $\mathcal{S}$. From Lemma 6 and we use the inequality
$\sigma_{c^{\prime}}({\bm{B}}_{0})>\|\hat{\bm{\Delta}}\|_{2}$ which ensures
full-rankness of $\hat{{\bm{B}}}^{\ast}$, which is the key ingredient in order
to prove that
$\operatorname*{span}(\hat{{\bm{B}}}^{\ast})=\mathcal{S}_{\perp}$. We can then
Use Theorem 9 for matrix ${\bm{B}}_{0}$. Note that columns of ${\bm{B}}_{0}$
are drawn independently and are uniformly distributed on the unit sphere.
Hence, columns of ${\bm{B}}_{0}$ are sampled by subgaussian distribution and
the LHS of the inequality of the theorem appears if we scale with
$\frac{1}{\sqrt{D}}$ so that to create unit-norm columns and use LHS of the
inequality of Theorem 7. The RHS of the inequality is due to the upper bound
of $\|\hat{\bm{\Delta}}\|_{2}$ as stated in Lemma 8. The absolute constants
$C_{1},C_{2}$ depend only the subgaussian norm of the uniform distribution
(they is no dependency on the dimensions of the problem).

##  Appendix B Experimental details and additional material

All experiments were conducted on a MacBook Pro 2.6GhZ 6-Core Intel Core i7,
memory 16GB 2667 Mhz DDR using Matlab2019B. For computational purposes and in
order to avoid fine-tuning of the piecewise geometrically diminishing (PGD)
step size, the modified backtracking line-search (MBLS) step-size rule was
adopted for DPCP-PSGM as proposed in Zhu et al. (2018). We define the distance
between two subspaces spanned by matrices ${\bm{B}}$ and ${\bm{A}}$ as
$dist({\bm{B}},{\bm{A}})=\min_{{\bm{Q}}\in\mathbb{O}(D,c)}\|{\bm{B}}-{\bm{A}}{\bm{Q}}\|_{F}$
where $\mathbb{O}(D,c)$ denotes the Stiefel manifold of orthogonal matrices of
rank $c$. Note that
$dist({\bm{B}},{\bm{A}})=0\iff\operatorname*{span}({\bm{B}})\equiv\operatorname*{span}({\bm{A}})$
(see Zhu et al. (2019)).

###  B.1 Outliers Pursuit in Washington DC Mall AVIRIS HSI

Hyperspectral images (HSIs) provide rich spectral information as compared to
RGB images capturing a wide range of the electromagnetic spectrum. Washington
DC Mall AVIRIS HSI contains contiguous spectral bands captured at 0.4 to
2.4$\mu$m region of visible and infrared spectrum, Giampouras et al. (2019).
In this experiments we randomly choose 10 out of its 210 spectral bands. Due
to high coherence in the both the spectral and the spatial domain, pixels of
HSIs admit representations in low-dimensional subspaces. Here, we use a
100$\times$100 segment of the hyperspectral image selecting randomly 10 out of
its $D=210$ spectral bands. We form a matrix $\bm{\tilde{\mathcal{X}}}$ of
size $10\times 10000$ whose columns correspond to different points in the
$10$-dimensional ambient space. Then we corrupt columns of
$\bm{\tilde{\mathcal{X}}}$ by replacing them with outliers that are generated
uniformly at random with unit $\ell_{2}$ norm for two different outliers’
ratios i.e., $r=0.8$ and $r=0.9$. In the corrupted $\bm{\tilde{\mathcal{X}}}$,
the remaining clear pixels are considered as the inliers. Table 1 displays the
F1 scores obtained by DPCP-PSGM, RSGM and DPCP-IRLS algorithm. The latter two
algorithms are evaluated in two scenarios: a) codimension is initialized
$c^{\prime}=5$ and b) $c^{\prime}=10$. Given the singular value distribution
of the initial image, we infer that the dimension $d$ of the inliers’ subspace
is less or equal than 5.

Table 1: Results on Washinghton DC AVIRIS hyperspectral image Methods | F1-scores
---|---
| r = 80% | r = 90%
DPCP-PSGM (unknown $c$) | 0.994 | 0.993
RSGM (unknown $c$) | 0 | 0
DPCP-IRLS (unknown $c$) | 0 | 0
RSGM ($c=5$) | 0.999 | 0.993
DPCP-IRLS $c=5$ | 1 | 0.995

Hence, $c^{\prime}=5$ (recall $c=D-d$) is close to the true codimension value
while $c^{\prime}=10$ is an overestimate thereof. From Table 1, we can see
that the proposed DPCP-PSGM succeeds in both outliers’ ratios regardless its
unawareness of the true codimension value. On the other hand, DPCP-IRLS and
RSGM fail when initialized with $c=10$ and this is attributed to the
restrictions induced due to the orthogonality constraints they both impose. In
Fig. 4 we provide annotated versions of the clean HSI, its corrupted by
outliers version for outliers’ ratio $r=90\%$, and the annotated outliers as
recovered by the proposed DPCP-PSGM, RSGM, RSGM with $c^{\prime}=5$ and DPCP-
IRLS with $c^{\prime}=5$.

|  |
---|---|---
(a) | (b) | (c)
|  |
(d) | (e) | (f)
Figure 4: (a) False RGB color image of the clean version of Washington Mall
AVIRIS HSI, (b) corrupted by outliers depicted with red and inliers
correpsonding the non-red pixels (c) annotated outliers as recoverd by the
proposed DPCP-PSGM method initialized with $c^{\prime}=10$ (d) RSGM with
$c^{\prime}=10$, (e) RSGM with $c^{\prime}=5$ and (f) DPCP-IRLS with
$c^{\prime}=5$.