# Provable Low Rank Plus Sparse Matrix Separation Via Nonconvex
Regularizers††thanks: This work was supported in part by National Science
Foundation under Grant Number DMS-1736326

April Sagan Department of Biomedical Informatics, School of Medicine,
University of Pittsburgh, Pittsburgh, USA; UPMC Hillman Cancer Center,
University of Pittsburgh, Pittsburgh, USA (,http://www.aprilsagan.net). John
E. Mitchell Department of Mathematical Sciences, Rensselaer Polytechnic
Institute, Troy, NY ()

###### Abstract

This paper considers a large class of problems where we seek to recover a low
rank matrix and/or sparse vector from some set of measurements. While methods
based on convex relaxations suffer from a (possibly large) estimator bias, and
other nonconvex methods require the rank or sparsity to be known a priori, we
use nonconvex regularizers to minimize the rank and $l_{0}$ norm without the
estimator bias from the convex relaxation. We present a novel analysis of the
alternating proximal gradient descent algorithm applied to such problems, and
bound the error between the iterates and the ground truth sparse and low rank
matrices. The algorithm and error bound can be applied to sparse optimization,
matrix completion, and robust principal component analysis as special cases of
our results.

Keywords: Nonconvex Regularizers, Low Rank Models, Sparse Optimization, Matrix
Completion, Robust PCA

AMS subject classifications: 62J07, 15A83, 90C26

##  1 Introduction

In order to better understand large datasets and to make inferences about
them, it is helpful to understand the underlying patterns in the datasets.
Even when the underlying pattern is highly nonlinear, the data matrix can be
approximated as being low rank, an observation that enables techniques to
analyze the data in terms of a low dimensional latent space, such as Principal
Component Analysis (PCA), identifying outliers through Robust PCA (RPCA), and
accurately inferring data points from very few observations of a data matrix
through matrix completion.

Data analysis techniques based upon this low rank property have received much
attention in the past decade, with impressive computational results on large
matrices and theoretical results guaranteeing the success of RPCA and matrix
completion [8][4]. Many of these results are based on minimizing the nuclear
norm of a matrix (defined as the sum of the singular values) as a surrogate
for the rank function, similar to minimizing the $l_{1}$ norm to promote
sparsity in a vector.

While the convex relaxation is an incredibly useful technique in many
applications, minimizing the nuclear norm of a matrix has been shown to
introduce a (sometimes very large) estimator bias. Intuitively, we expect to
see this bias because if we hope to recover a rank $r$ matrix, we must impose
enough weight on the nuclear norm term so that the $(r+1)$th singular value is
zero. By the nature of the nuclear norm, this requires also putting weight on
minimizing the first $r$ singular values, resulting in a bias towards zero
proportional to the spectral norm of the noise added to the true data matrix.

Fortunately, recent work has shown that the estimator bias from convex
regularizers can be reduced (or even eliminated, for well conditioned
matrices) by using nonconvex regularizers, such as the Schatten-p norm or the
minimax concave penalty (MCP). It has been shown that for sparse optimization,
the nonconvexity introduced from these regularizers does not create a further
burden in the the optimization process – in the right circumstances, the
nonconvex problem has just one minimizer [21]. Similar results for rank
minimization problems have been previously unavailable, a gap that we have
aimed to fill in this paper.

###  1.1 Summary of Contributions

In this paper, we focus on the nonconvex, unconstrained optimization problem
where we find a low-rank matrix $L\in\mathbb{R}^{d_{1}\times d_{2}}$ and
sparse vector $s\in\mathbb{R}^{d_{s}}$.

$\underset{L,s}{\text{min
}}\;\frac{\lambda_{L}}{d_{1}d_{2}}\Phi_{\gamma_{L}}(L)+\frac{\lambda_{s}}{d_{s}}\phi_{\gamma_{s}}(s)+\frac{1}{2n}||\mathcal{A}_{L}(L)+A_{S}s-b||_{2}^{2}$
(1)

The linear mappings $\mathcal{A}_{L}:\mathbb{R}^{d_{1}\times
d_{2}}\rightarrow\mathbb{R}^{n}$ and $A_{s}\in\mathbb{R}^{d_{s}\times n}$
serve as the observation models of the underlying low rank matrices and sparse
vectors. Most commonly, we are interested in the observation model

$\big{[}\mathcal{A}_{\Omega^{obs}}(X)\big{]}_{k}=X_{i_{k},j_{k}}$

for $(i_{k},j_{k})\in\Omega^{obs}$, where $\Omega^{obs}\subseteq\\{1,\ldots
d_{1}\\}\times\\{1,\ldots d_{2}\\}$ is the set of indices where we have a
measurement of the low rank matrix we hope to reconstruct.

We denote $\phi:\mathbb{R}\rightarrow\mathbb{R}_{+}$ to be a concave function
used to promote sparsity in both the singular values of $L$ and individual
entries in $s$. We overload the notation to allow for $\phi_{\gamma_{s}}$ to
be a function of a vector $x\in\mathbb{R}^{d_{s}}$ whose range is
$\mathbb{R}_{+}$, and we denote $\Phi_{\gamma_{L}}:\mathbb{R}^{d_{1}\times
d_{2}}\rightarrow\mathbb{R}_{+}$ as a surrogate to the rank function:

$\phi_{\gamma_{s}}(x)=\sum_{i=1}^{d_{s}}\phi_{\gamma_{s}}(x_{i}),\;\>\Phi_{\gamma_{L}}(X)=\sum_{i=1}^{\text{min}(d_{1},d_{2})}\phi_{\gamma_{L}}(\sigma_{i}(X)).$

where $\sigma_{i}(X)$ denotes the $i$th largest singular value of $X$. We
restrict our focus to nonconvex regularizers that are amenable regularizers,
as described in [20], and defined below.

######  Definition 1.

A function $\phi_{\gamma}(t):\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is
amenable if it satisfies the following criteria.

  1. 1.

$\phi_{\gamma}(0)=0$

  2. 2.

$\phi_{\gamma}$ is non decreasing

  3. 3.

For $t>0$, the function $\frac{\phi_{\gamma}(t)}{t}$ is non increasing in $t$.

  4. 4.

the function $\phi_{\gamma}$ is differentiable for all $t\neq 0$ and
subdifferentiable at $t=0$ with $\lim_{t\rightarrow
0^{+}}\phi_{\gamma}^{\prime}(t)=1$.

  5. 5.

The function $\phi_{\gamma}(t)$ is $\nu$ weakly convex. That is, the function
$\rho_{\nu}:=\phi_{\gamma}(t)+\frac{\nu}{2}t^{2}$ is convex.

We present a very simple alternating algorithm to find a stationary point of
(1). Though similar algorithms have been previously studied and shown to
converge, we present a novel analysis of this algorithm and show that not only
does this algorithm linearly converge, but it converges to the exact low rank
matrix $L^{*}$ and sparse vector $s^{*}$ when no Gaussian noise is present.
Furthermore, when Gaussian noise is present, we obtain an error bound that
matches the minmax optimal rate. Our bound greatly improves upon bounds
obtained in previous results analysing the convex relaxation and quantify the
observation based on computational results that nonconvex regularizers reduce
the impact of noise on the quality of the estimator.

###  1.2 Related Works

The matrix completion problem is as follows: given the values of a low rank
matrix for only a sparse set of indices, we seek to determine the rest of the
values of the matrix. While the problem of finding the minimum rank matrix
that fits the observations is NP hard in general, it has been shown that under
some assumptions, the global minimizer to the convex problem

$\underset{X\in\mathbb{R}^{d_{1}\times
d_{2}}}{\text{min}}||X||_{*}\;\text{s.t.}\;X_{ij}=M_{ij}\;\forall(i,j)\in\Omega^{obs}$

is exactly $M$, where $\Omega^{obs}$ is the set of indices of $M$ we have
observed. If $M$ is rank $r$ and $\mu-$incoherent (as defined in section 3.1),
then with high probability for a set, $\Omega^{obs}$ of $n$ indices chosen
uniformly at random, $M$ is the unique minimizer to the convex relaxation so
long as $n>C\mu rd_{1}\log^{1.2}(d_{1})$ for some universal constant $C$ [8,
4]. This condition was later improved to $n>C\mu rd_{1}\log(d_{1})$ [29].

Using a convex relaxation for Robust PCA has similar results. While PCA is a
powerful technique, it has been shown to be less reliable when just a sparse
set of data points are grossly corrupted, and so the goal of RPCA is to
identify and remove such corruptions by separating the data matrix into the
sum of a low rank and sparse matrix.

$\underset{L,s}{\text{min}}{||L||_{*}+\lambda_{0}||S||_{1}}\>\text{ st
}\>L_{ij}+S_{ij}=M_{ij}\;\forall(i,j)\in\Omega^{obs}$

The convex relaxation was shown to give the exact solution when every entry of
$M$ is observed in [10], and when only partially observed (under the same
assumptions necessary in matrix completion) by [6]. In contrast to this paper,
these works assume that there is measurement noise (besides the sparse
corrupted entries).

In many cases of practical interest, our measurements may have some level of
noise in addition to being only partially observed or having some corrupted
entries. For matrix completion, we can relax the constraint using a penalty
formulation as follows:

$\underset{X\in\mathbb{R}^{d_{1}\times
d_{2}}}{\text{min}}\lambda||X||_{*}+\sum_{(i,j)\in\Omega^{obs}}\big{(}X_{ij}-M_{ij}\big{)}^{2}$

Likewise, RPCA can be formulated as solving:

$\displaystyle\underset{L,S\in\mathbb{R}^{d_{1}\times d_{2}}}{\text{min}}$
$\displaystyle{\lambda_{L}||L||_{*}+\lambda_{s}||S||_{1}}+\sum_{(i,j)\in\Omega^{obs}}\big{(}L_{ij}+S_{ij}-M_{ij}\big{)}^{2}$
(2)

Statistical guarantees on the performance of the first of these estimators are
discussed in [7, 26, 25], and the latter in [2]. Specific bounds and
assumptions for these works are discussed in Section 4.

In order to reduce the estimator bias for $l_{1}$ minimization, [9] proposed
an iteratively reweighted $l_{1}$ norm method to place more weight on
minimizing smaller entries, and less on entries further from zero. This idea
was generalized to minimizing any amenable regularizes to promote sparsity.
Theoretical results on the subject include algorithmic guarantees similar to
the ones presented in this paper [34], and a proof that the nonconvex problem
has no spurious local minimizers [20, 21].

The same regularizers could be used as a surrogate to the rank function, as
originally proposed by [24]. In [22, 23], the authors propose a generalization
of the singular value thresholding algorithm proposed by [3], which was later
applied to the problem of RPCA in [11, 18]. For the problem of matrix
completion, the algorithms proposed by [36] and [30] achieve the fastest
computational complexity in other state of the art methods.

Other approaches to low rank optimization rely on, instead of minimizing a
surrogate to the rank function, constraining the matrix to be a given rank.
This can be done using constrained optimization to optimize over the set of
all rank $r$ matrices [33], or by using the low-rank factorization of a matrix
$X=UV^{T}$ for $U\in\mathbb{R}^{d_{1}\times r}$ and
$V\in\mathbb{R}^{d_{2}\times r}$ [35, 31].

Previous work on theoretical results pertaining to rank and sparsity
constrained methods have consisted of algorithmic guarantees ensuring that we
can obtain a matrix sufficiently close to the ground truth low rank and/or
sparse matrix for both matrix completion [17] and RPCA [27, 38, 37].
Additionally, both of these problems have been shown to have no spurious local
minimizers, and so the ground truth matrices are the only minimizers under
some assumptions [13] [12].

Table 1: Table of common nonconvex regularizers used for $l_{0}$ and rank minimization, along with the associated proximal operator. | $\phi_{\gamma}(x)$ | $\text{prox}_{\phi_{\gamma}}^{\tau}(y)$
---|---|---
$l_{1}$ | $|x|$ | $\text{sign}(y)(|y|-\frac{1}{\tau})_{+}$
| Capped
---
$l_{1}$ Norm
$\text{min}(|x|,\frac{c}{2})_{+}$ | $\begin{cases}\text{sign}(y)(|y|-\frac{1}{\tau})&|y|\leq\frac{c}{2}\\\ y&\text{otherwise }\end{cases}$
SCAD | $\begin{cases}|x|&|x|\leq 1\\\ \frac{2\alpha|x|-x^{2}-1}{2(\alpha-1)}&1\leq|x|\leq\alpha\\\ \frac{(\alpha+1)}{2}&|x|>\alpha\end{cases}$ |  $\frac{1}{\alpha-1}$ $\begin{cases}\text{sign}(y)(|y|-\frac{1}{\tau})_{+}&|y|\leq 1\\\ \frac{\text{sign}(y)\tau(a-1)}{\tau(a-1)-1}(|y|-\frac{a}{\tau})&\leq|y|\leq\alpha\\\ y&|y|>\alpha\end{cases}$
MCP | $\begin{cases}|x|-\frac{x^{2}}{2}&|x|\leq 1\\\ \frac{1}{2}&|x|>1\end{cases}$ | $\begin{array}[]{l}\text{sign}(y)\text{min}\bigg{(}|y|,\hskip 14.22636pt\frac{\tau}{\tau-1}(|y|-\frac{1}{\tau})_{+}\bigg{)}\end{array}$

##  2 Alternating Proximal Gradient Descent Algorithm

Many different methods have been shown to be effective when minimizing
nonconvex relaxations of the $l_{0}$ and rank functions, including iterative
reweighted methods, and methods based on low rank factorization. In this
paper, we focus on the most commonly used technique: alternating proximal
gradient descent.

Consider the objective function $F(x)=\phi(x)+g(x)$, where $g(x)$ is convex
and differentiable, and $\phi(x)$ is weakly convex. Instead of minimizing
$F(x)$ directly, the proximal gradient descent method approximates $g(x)$ by a
quadratic, strongly convex function $\bar{g}_{k}(x)$ centered about the point
$x^{k}$.

$\bar{g}_{k}(x)=g(x^{k})+\langle\nabla
g(x^{k}),x-x^{k}\rangle+\frac{1}{2\tau}||x-x^{k}||^{2}$

At each iteration, we now minimize the function $F_{k}(x)=\phi(x)+g_{k}(x)$.
For sufficiently small $\tau$, this function is strongly convex, and the
proximal gradient descent algorithm is guaranteed to converge.

The proximal gradient descent method applied to the function
$F(x)=\phi(x)+\bar{g(x)}$ iteratively solves the following problem:

$\displaystyle x^{k+1}=\underset{x}{\text{argmin
}}\phi(x)+\frac{1}{2\tau}||x-\big{(}x^{k}-\tau\nabla
g(x^{k})\big{)}||^{2}:=\text{prox}_{\phi}^{\tau}\big{(}x^{k}-\tau\nabla
g(x^{k})\big{)}$

where we defined the proximal operator of a function as the minimum of a
combination of the function and the distance from a given point. For many
functions $\phi$ that we are interested in, the proximal operator has a closed
form solution, some of which are shown in Table 1.

For each of the sparsity promoting regularizers in Table 1, the proximal
operator is also dubbed as a shrinkage operator or a thresholding operator
because when the input is less than $\frac{1}{\tau}$, the output is 0.
Otherwise, the input is moved towards zero or, for some nonconvex
regularizers, is unchanged. So, we can view the proximal gradient algorithm as
iteratively taking a step in the gradient direction of $g(x)$, and then
applying the proximal operator to promote sparsity.

When applied to the optimization problem in Equation (1), we have

$\displaystyle\tilde{L}^{k+1}=$ $\displaystyle
L^{k}-\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\big{(}\mathcal{A}_{L}(L^{k})+A_{s}s^{k}-b\big{)}$
(3a) $\displaystyle L^{k+1}=$
$\displaystyle\text{prox}_{\Phi_{\gamma_{L}}}^{{\tau_{L}}{\lambda_{L}}}(\tilde{L}^{k+1})$
(3b) $\displaystyle\tilde{s}^{k+1}=$ $\displaystyle
s^{k}-{\tau_{S}}\frac{d_{s}}{n}A_{s}^{T}\big{(}\mathcal{A}_{L}(L^{k+1})+A_{s}s^{k}-b\big{)}$
(3c) $\displaystyle s^{k+1}=$
$\displaystyle\text{prox}_{\phi_{\gamma_{s}}}^{{\tau_{s}}{\lambda_{s}}}(\tilde{s}^{k+1})$
(3d)

To further simplify the problem, the following proposition will allow $L$
subproblem to be solved in each singular value separately.

######  Proposition 1.

Consider the optimization problem

$\underset{X}{\text{min
}}\sum_{i}\phi_{\gamma}(\sigma_{i}(X))+\frac{1}{2\tau}||X-Y||_{F}^{2}$ (4)

where $\phi_{\gamma}$ is a $\nu$ weakly convex function. If $\nu<\tau$, then
equation (4) is strongly convex and the minimizer $X^{*}$ has the same
singular vectors as $Y$ with singular values given by

$\displaystyle\sigma_{i}(X^{*})=$ $\displaystyle\underset{x}{\text{argmin
}}\phi_{\gamma}(x)+\frac{1}{2\tau}(x-\sigma_{i}(Y))^{2}$ $\displaystyle:=$
$\displaystyle\text{prox}_{\phi_{\gamma}}^{\tau}(\sigma_{i}(Y))$

where $\text{prox}_{\phi_{\gamma}}^{\tau}(\sigma_{i}(Y))$ is the proximal
operator.

###### Proof.

By convexity of $\phi$,

$\phi(\lambda_{i}(X_{2}))\geq\phi(\lambda_{i}(X_{1}))+\phi^{\prime}(\lambda_{i}(X_{1}))(\lambda_{i}(X_{2})-\lambda_{i}(X_{1}))$

Summing over all $i=1,\ldots n$,

$\Phi(X_{2})\geq\Phi(X_{1})+\langle\nabla\phi(\lambda(X_{1})),\lambda(X_{1})-\lambda(X_{2})\rangle$

By Corollary 1,

$\Phi(X_{2})\geq\Phi(X_{1})+\langle
U_{1}\nabla\phi(\lambda(X_{1}))U_{1}^{T},X_{1}-X_{2}\rangle$

And, as $\nabla\Phi(X)=U\nabla\phi(\lambda(X))U^{T}$, $\Phi$ is convex.

Consider if the global optimizer $X^{*}$ was not of the form $U\Sigma V^{T}$,
for a diagonal matrix $\Sigma$. Let $\bar{X}=U\Sigma^{*}V^{T}$, where
$\Sigma^{*}_{ii}=\sigma_{i}(X^{*})$. By the of Hoffman-Wielandt inequality
(Corollary A.1),

$||X^{*}-Y||^{2}_{F}\leq\sum_{i}\big{(}\sigma_{i}(X^{*})-\sigma_{i}(Y^{*})\big{)}^{2}=||X^{*}-Y||^{2}.$

And, because $\Phi(\bar{X})=\Phi(X^{*})$, $\bar{X}$ must also be a global
minimizer, contradicting the premise that $X^{*}$ is the sole global
minimizer. ∎

Proposition 1 tells us that $L^{k}$ has the singular vectors of
$\tilde{L}^{k}$ and singular values given by

$\sigma_{i}(L^{k})=\text{prox}_{\phi_{\gamma}}^{{\tau_{L}}{\lambda_{L}}}(\sigma_{i}(\tilde{L}^{k})).$

Likewise, the subproblem in $S$ can be solved in each entry of $s$
individually.

Algorithm 1 Alternating Proximal Gradient Descent for Low-Rank Plus Sparse
Optimization (APGD)

for $k=1,\ldots$ do

$g^{k+1}=-\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\big{(}\mathcal{A}_{L}(U^{k}\Sigma^{k}V^{k})+A_{s}s^{k}-b\big{)}$

$[U^{k+1},\tilde{\Sigma}^{k+1},V^{k+1}]=\text{LRSSVD}(U^{k},{\Sigma}^{k},V^{k},g^{k+1})$

${\Sigma}^{k+1}=\text{prox}\big{(}\tilde{\Sigma}^{k+1}\big{)}$

$\tilde{s}^{k}=s^{k}-\tau_{s}\frac{d_{s}}{n}A_{s}^{T}\big{(}\mathcal{A}_{L}(U\Sigma
V)+A_{s}s^{k}-b\big{)}$

$s^{k+1}=\text{prox}^{\lambda_{s}}\big{(}\tilde{s}^{k+1}\big{)}$

end for

The slowest operation in the alternating proximal gradient method is by far
the singular value decomposition. However, in practice we can reduce the
number of operations by calculating the truncated singular value decomposition
only using the first $r_{0}$ singular values, where $r_{0}$ is an upper bound
on the rank, and enforce that the remaining singular values are zero.
Alternatively, we can calculate each singular value in descending order and
stop when a singular value falls below $\lambda_{L}$, as all remaining
singular values will be set to zero by the proximal operator. So, in the case
of RPCA where each entry is observed, each iteration has a computational
complexity of $\mathcal{O}(d_{1}d_{2}r_{0})$, which matches other state of the
art methods.

In the case of matrix completion, however, only a sparse set of entries of the
low rank matrix are observed, which could be used to increase the efficiency
by reducing the amount of computation needed to find the singular value
decomposition of $L$ at each iterations. For a low rank matrix with a low rank
factorization $U\Sigma V$, we refer to the problem of finding the SVD of the
rank $r$ approximation to the matrix $U\Sigma V+g$ for a sparse matrix $g$ as
Low Rank plus Sparse SVD (LRSSVD), originally proposed by [16].

The LRSSVD task can be accomplished efficiently using the same methods as if
we were to find the SVD of any other matrix, such as the Power Iteration
method. Recall that the computational complexity of the Power Iteration is
limited by the amount of operations needed to multiply the matrix by a vector.
Because the computational complexities of calculating both $u(U\Sigma
V^{T}+Y)$ for $u\in\mathbb{R}^{d_{1}}$ and $(U\Sigma V^{T}+Y)v$ for
$v\in\mathbb{R}^{d_{2}}$ are $\mathcal{O}((d_{1}+d_{2})r_{0}+n)$, we can
calculate the top $r_{0}$ singular values and vectors of $X+Y$ with only
$\mathcal{O}\big{(}(d_{1}+d_{2})r_{0}^{2}+nr_{0}\big{)}$ operations.

Algorithm 2 Singular Value Decomposition for a Low Rank Plus Sparse Matrix

0: $U,\Sigma,V,Y$

0: Singular value decomposition of $U\Sigma V^{T}+Y$ Initialization:
$\tilde{V}=V$

1: for $k=1,..,$ do

2: $\tilde{U}^{k}=(U\Sigma
V^{T}+Y)\tilde{V}(\tilde{V}\tilde{V}^{T})^{\dagger}$

3: $\tilde{V}^{k}=(\tilde{U}\tilde{U}^{T})Z^{\dagger}\tilde{U}^{T}(U\Sigma
V^{T}+Y)$

4: end for

5: $[Q^{U},R^{U}]=\text{QR}(\tilde{U}^{k})$,
$[Q^{V},R^{V}]=\text{QR}(\tilde{V}^{k})$

6: $[U^{R},\Sigma^{R},V^{R}]=\text{SVD}(R^{U}\tilde{\Sigma}{R^{V}}^{T})$

7: $\tilde{U}=Q^{U}U^{R}$, $\tilde{V}=Q^{V}V^{R}$Return:
$(\tilde{U},\tilde{\Sigma},\tilde{V})$

The other operations in Algorithm 1 take no more time than the LRSSVD. The
gradient in the $L$ direction, $g^{k}$, requires calculating
$\Sigma_{ii}^{k}U^{k}_{i}(V_{i}^{k})^{T}$ for each entry in the support of
$\mathcal{A}_{L}$. In the case of matrix completion, this is $nr$ operations,
which matches the computational complexity per iteration for state of the art
matrix completion algorithms.

##  3 Analysis of APGD Algorithm

In this section, we present the main result of the paper: a recursive bound on
the difference of the iterates of the alternating proximal gradient algorithm
and the ground truth low rank matrix and sparse vector. We present the bound
for the most general case, and give results on specific problems in the
following section.

###  3.1 Restricted Isometry and Orthogonality Properties

In order to bound the error in the output of our algorithm relative to the
underlying ground truth low-rank and sparse matrices $L^{*}$ and $s^{*}$, we
must first make a number of assumptions about $L^{*}$, $s^{*}$, and the
observation models $\mathcal{A}_{L}$ and $A_{s}$.

First, we must assure that the low rank matrix $L^{*}$ can be separated from a
sparse matrix – that is, $L^{*}$ is not sparse itself. Not only is this
necessary for low-rank plus sparse decomposition, but for the problem of
matrix completion, this assumption is necessary to assure that a sparse set of
observations is a good representation of the entire matrix. For example,
consider the matrix consisting of zeros in every entry besides one entry,
where the value is 1. We must observe every entry in the matrix to assure that
we can reconstruct the matrix exactly, due to the fact that we must observe
the nonzero entry and every entry in its row and column. To exclude such ill-
posed problems from our analysis, we will assume that $L^{*}$ is incoherent,
defined as:

######  Definition 2 ([8]).

Let $X\in\mathbb{R}^{d_{1}\times d_{2}}$ be a rank $r$ matrix with singular
value decomposition $X=U\Sigma V$, for orthonormal matrices
$U\in\mathbb{R}^{d_{1}\times r},V\in\mathbb{R}^{d_{2}\times r}$, and diagonal
matrix $\Sigma\in\mathbb{R}^{r\times r}$. The tangent space of $X$ is defined
as

$\mathcal{T}=\bigg{\\{}UA+BV\bigg{|}A\in\mathbb{R}^{r\times
d_{2}},B\in\mathbb{R}^{d_{1}\times r}\bigg{\\}}.$ (5)

Furthermore, we say the matrix $X$ (or its tangent space $\mathcal{T}$) is
$\mu$ incoherent if

$||U_{:i}||^{2}_{2}\leq\frac{\mu r}{d_{1}},\;||V_{j:}||^{2}_{2}\leq\frac{\mu
r}{d_{2}}\;\;\;\forall i,j$ (6)

We define the projection of a matrix onto the sparse space $\Omega$ as

$\mathcal{P}_{\Omega}\big{(}X\big{)}=\begin{cases}X_{ij}&(i,j)\in\Omega\\\
0&\text{else}\end{cases}$ (7)

and onto the tangent space $\mathcal{T}$ as

$\mathcal{P}_{\mathcal{T}}\big{(}X\big{)}=UU^{T}X+XVV^{T}-UU^{T}XVV^{T}$ (8)

Next, we discuss the conditions that the observation models $\mathcal{A}_{L}$
and $A_{s}$ must satisfy in order to recover the ground truth low rank and
sparse matrix, known as the restricted isometry property. Loosely, the RIP
states that for any two vectors in $\Omega$ (or matrices in $\mathcal{T}$), we
can obtain a sufficiently accurate estimate of the distance between the two
through the observation model $A_{s}$ (or $\mathcal{A}_{L}$).

The RIP was originally proposed for sparse vectors by Candès and Tao in [5].
Original attempts to extend the property to low rank matrices failed to
consider coherent matrices, and thus had little practical applications. Candès
and Tao [8] later introduced the incoherence assumption and proved that it
applied to the problem of matrix completion. Here, we give two versions of a
definition of restricted isometry property, one for sparse vectors and one for
matrices that have a low-rank, incoherent tangent space.

######  Definition 3.

The linear mapping $A_{s}$ satisfies the $(\alpha,\kappa)$ sparse Restricted
Isometry Property if, for any $x$ satisfying $||x||_{0}\leq\alpha d_{s}$

$\displaystyle(1-\kappa_{S})||x||^{2}$
$\displaystyle\leq\tau_{s}\frac{n}{d_{s}}||A_{s}x||^{2}\leq(1+\kappa_{S})||x||^{2}$

for some constant $\tau_{s}$. Likewise, the linear mapping $\mathcal{A}$
satisfies the $(\mu,r,\kappa)$ low rank Restricted Isometry Property if for
any $X$ in a $\mu$-incoherent rank $r$ tangent space $\mathcal{T}$,

$\displaystyle(1-\kappa_{L})||X||_{F}^{2}$
$\displaystyle\leq\tau_{L}\frac{n}{d_{1}d_{2}}||\mathcal{A}_{L}X||_{F}^{2}\leq(1+\kappa_{L})||X||_{F}^{2}$

for some constant $\tau_{L}$.

In some cases, it may be more useful to use the following characterization of
the RIP, which bounds the difference between the operator
$\tau\mathcal{A}^{*}\mathcal{A}$ and the identity operator when restricted to
sparse vectors or low-rank and incoherent matrices.

######  Proposition 2.

For a matrix $A_{s}\in\mathbb{R}^{n\times d_{s}}$ satisfying the
$(\alpha,\kappa_{s})$ sparse RIP,

$\displaystyle||\frac{\tau_{s}n}{d_{s}}\mathcal{P}_{\Omega}A_{s}^{T}A_{s}\mathcal{P}_{\Omega}-\mathcal{P}_{\Omega}||^{2}\leq\kappa_{s}.$

Likewise, for any linear mapping $\mathcal{A}_{L}:\mathbb{R}^{d_{1}\times
d_{2}}\rightarrow\mathbb{R}^{n}$ satisfying the $(\mu,r,\kappa_{L})$ low rank
RIP,

$\displaystyle||\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}\mathcal{A}_{L}\mathcal{P}_{\mathcal{T}}-\mathcal{P}_{\mathcal{T}}||^{2}\leq\kappa_{L}$

Finally, we discuss the interplay between the set of sparse matrices and the
low rank, incoherent tangent space, and their observation models. We hope to
be able to separate the measurement vector $b$ into two parts: one in the span
of $\mathcal{A}_{L}\mathcal{P}_{\mathcal{T}}$, and one in the span of
$A_{s}\mathcal{P}_{\Omega}$. In order to achieve this quickly, we require that
there are no non-trivial vectors in the intersection of the two sets, which is
equivalent to saying that
$||\mathcal{P}_{\Omega}A_{s}^{T}\mathcal{A}_{L}\mathcal{P}_{\mathcal{T}}||<1$.
Under some assumptions, this norm is actually close to zero, a concept we
refer to as restricted orthogonality, which we define here and verify that it
applies to the problems we are interested in in Section 4.

######  Definition 4.

Linear maps $\mathcal{A}_{L}$ and $A_{s}$ satisfy the $\kappa-$ restricted
orthogonality property over the sets $\mathcal{T}$ and $\Omega$ (respectively)
when

$\displaystyle\frac{d_{s}}{n}||\mathcal{P}_{\Omega}A_{s}^{*}\mathcal{A}_{L}\mathcal{P}_{\mathcal{T}}||^{2}\leq\kappa,\;\frac{d_{1}d_{2}}{n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}A_{s}\mathcal{P}_{\Omega}||^{2}\leq\kappa$

###  3.2 Main Result

Define the difference between the iterates of the alternating proximal
gradient descent algorithm and the ground truth low rank matrix and sparse
vector at iteration $k$ as $\Delta_{L}^{k}=L^{*}-L^{k}$ and
$\Delta_{s}^{k}=s^{*}-s^{k}$. Our main result in the most general form gives a
bound on the norm of $\Delta_{L}^{k}$ and $\Delta_{s}^{k}$ in terms of the
differences at the previous iteration, $\Delta_{L}^{k-1}$ and
$\Delta_{s}^{k-1}$.

######  Theorem 1.

Let $L^{k}$ and $s^{k}$ be the sequences generated by Algorithm 1. Assume that

$b=\mathcal{A}_{L}(L^{*})+A_{s}s^{*}+\mathcal{E}\in\mathbb{R}^{n},$

where $L^{*}\in\mathbb{R}^{d_{1}\times d_{2}}$ is a rank $r$ and $\mu-$
incoherent matrix, and $s^{*}\in\mathbb{R}^{d_{s}}$ is a sparse vector with
$supp(s^{*})=\Omega$, and the linear mappings $\mathcal{A}_{L}$ and $A_{s}$
satisfy the $(2r,3\mu,\kappa_{L})-$ low rank RIP and the
$(\alpha,\kappa_{s})-$ sparse RIP respectively, and together satisfy the ROP
with constant $\kappa$. If
$\lambda_{L}\geq||\mathcal{A}_{L}^{*}(\mathcal{E})||_{2}+||\mathcal{A}_{L}^{*}A_{s}\Delta_{s}^{k-1}||_{2}$,
then

$\displaystyle||\Delta_{L}^{k+1}$
$\displaystyle||_{F}^{2}\leq\kappa_{L}||\Delta_{L}^{k}||_{F}^{2}+\kappa\tau_{L}||\Delta_{s}^{k}||^{2}+\frac{\tau_{L}d_{1}d_{2}}{n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}\mathcal{E}||_{F}^{2}+\lambda_{L}r\phi^{\prime}(\sigma_{r}(L^{*}))$

Likewise, if
$\lambda_{S}\geq||A_{s}^{T}(\mathcal{E})||_{\infty}+||A_{s}^{T}\mathcal{A}_{L}\Delta_{L}^{k-1}||_{\infty}$
and $s_{min}$ is the smallest non-zero value of $s^{*}$, then

$\displaystyle||\Delta_{s}^{k+1}||^{2}$
$\displaystyle\leq\kappa_{s}||\Delta_{s}^{k}||^{2}+\kappa\tau_{s}||\Delta_{L}^{k+1}||_{F}^{2}+\frac{\tau_{s}d_{s}}{n}||\mathcal{P}_{\Omega}A_{s}^{T}\mathcal{E}||^{2}+\lambda_{s}\alpha
d_{s}\phi^{\prime}(s_{min}-\lambda_{s})$

For each of these bounds, we can think of the third term as the estimation
error introduced by the noise, and the fourth term as the approximation error,
which accounts for the bias in the regularizer proportional to the derivative
of the regularizer. Previous results for the nuclear norm and $l_{1}$ norm
give similar bounds, but make the concession that the approximation error is
the dominating term. Under some circumstances, that term is equal to zero in
our bound.

###  3.3 Proof of Main Result

We start by presenting the the following two lemmas regarding the proximal
operator for the low rank regularizers, which we prove in the following
section.

######  Lemma 1.

Let $L^{*}$ be a rank $r$, $\mu$-incoherent matrix whose singular vectors form
the tangent space $T$, and $\bar{L}$ be defined as

$\bar{L}=\underset{L\in\mathbb{R}^{d_{1}\times
d_{2}}}{\text{argmin}}\lambda\Phi(L)+||L-L^{*}+\delta||_{F}^{2}$ (9)

where $\Phi$ is an at most $\frac{1}{\lambda}$ weakly convex regularizer
satisfying Assumption 1, and $\delta\in\mathbb{R}^{d_{1}\times d_{2}}$
satisfies $||\delta||_{2}\leq\lambda$. Define $\Delta_{L}=\bar{L}-L^{*}$.
Then,

$||\Delta_{L}||^{2}_{F}\leq
2||\mathcal{P}_{\mathcal{T}}\big{(}\delta\big{)}||^{2}_{F}+\lambda
r\phi^{\prime}\big{(}\sigma_{r}(L^{*})\big{)}$ (10)

In order to utilize the RIP and ROP conditions, we need to verify that
$L^{k+1}$ is low rank and incoherent, and that $s^{k+1}$ is sparse, which we
will do inductively. Assume that $L^{k}$ is at most rank $r$, and that its
tangent space is $2\mu$ incoherent. Additionally, assume that
$\text{supp}\big{(}S^{k})\subseteq\text{supp}\big{(}S^{*})$. Clearly, these
conditions are met at the first iteration by initializing the algorithm with
$L^{0}=0$ and $S^{0}=0$.

At iteration $k+1$,

$\displaystyle L^{*}$ $\displaystyle-\bar{L}^{k+1}$ $\displaystyle=$
$\displaystyle
L^{*}-\big{(}L^{k}-\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\big{(}\mathcal{A}_{L}(L^{k})+A_{s}s^{k}-b)\big{)}\big{)}$
$\displaystyle=$ $\displaystyle
L^{*}-L^{k}+\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\mathcal{A}_{L}L^{k}+\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}A_{s}s^{k}-\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\big{(}\mathcal{A}_{L}(L^{*})+A_{s}(s^{*})+\mathcal{E}\big{)}$
$\displaystyle=$
$\displaystyle\big{(}\Delta_{L}^{k}-\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\mathcal{A}_{L}\Delta_{L}^{k}\big{)}-\tau_{L}\frac{d_{1}d_{2}}{n}\big{(}\mathcal{A}_{L}^{*}A_{s}\Delta_{S}^{k}-\mathcal{A}_{L}^{*}\mathcal{E}\big{)}$

where the first equality comes from the definition of $b$, and the second
inequality substitutes $\Delta_{S}^{k}=S^{*}-S^{k}$ and
$\Delta_{L}^{k}=L^{*}-L^{k}$. By Lemma 1 (along with the triangle inequality),
we have that

$\displaystyle||\Delta_{L}^{k+1}||_{F}^{2}\leq$
$\displaystyle||\mathcal{P}_{\mathcal{T}}\big{(}\Delta_{L}^{k}-\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\mathcal{A}_{L}\Delta_{L}^{k}\big{)}||_{F}^{2}+\tau_{L}\frac{d_{1}d_{2}}{n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}A_{s}\Delta_{S}^{k}||_{F}^{2}$
$\displaystyle+\tau_{L}\frac{d_{1}d_{2}}{n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}\mathcal{E}||_{F}^{2}+\lambda_{L}r\phi^{\prime}\big{(}\sigma_{r}(L^{*})-2{\lambda_{L}}{\tau_{L}}\big{)}$

Let $\mathcal{T}^{k}$ denote the union of the tangent space of the rank $r$
approximation $L^{k}$ and $\mathcal{T}$ so that
$\Delta_{L}^{k}\in\mathcal{T}^{k}$. Then,

$\displaystyle||\mathcal{P}_{\mathcal{T}}\big{(}\Delta_{L}^{k}-\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\mathcal{A}_{L}\Delta_{L}^{k}\big{)}||_{F}^{2}$
$\displaystyle\leq||\mathcal{P}_{\mathcal{T}^{k}}\big{(}\Delta_{L}^{k}-\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\mathcal{A}_{L}\Delta_{L}^{k}\big{)}||_{F}^{2}$
$\displaystyle\leq||\mathcal{P}_{T^{k}}\big{(}\tau_{L}\frac{d_{1}d_{2}}{n}\mathcal{A}_{L}^{*}\mathcal{A}_{L}-\mathcal{I}\big{)}\mathcal{P}_{\mathcal{T}^{k}}\Delta_{L}^{k}||_{F}^{2}\leq\kappa_{L}||\Delta_{L}^{k}||_{F}^{2}$

where the first inequality comes from the contractive property of
$\mathcal{P}_{\mathcal{T}}$, and the second inequality comes from the fact
that $\mathcal{P}_{\mathcal{T}^{k}}\Delta_{L}^{k}=\Delta_{L}^{k}$. And,
because $\mathcal{T}^{k}$ is has incoherence at most $3\mu$, we can apply the
low-rank RIP to obtain the third inequality.

By the inductive hypothesis stating that
$\text{supp}\big{(}S^{k})\subseteq\text{supp}\big{(}S^{*})$ (and thus, that
$S^{k}$ is $\alpha$-sparse), we can use the ROP to claim that

$\frac{d_{1}d_{2}}{n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}A_{s}\Delta_{S}^{k}||_{F}^{2}\leq\kappa||\Delta_{S}^{k}||_{F}^{2}.$

Combining these gives the desired bound on $\Delta^{k}_{L}$.

Now, we must show that $L^{k+1}$ is rank $r$ and $2\mu$ incoherent. By Weyl’s
inequality ([15], see Appendix A), we know that

$\displaystyle\sigma_{r+1}(\tilde{L}^{k+1})\leq$
$\displaystyle\sigma_{r+1}(L^{*})+||L^{*}-\tilde{L}^{k+1}||_{2}$
$\displaystyle\leq$
$\displaystyle||\Delta_{s}^{k}||_{2}+||\mathcal{A}_{L}^{*}\mathcal{E}||_{2}$

Because this is less than $\lambda_{L}$ by our assumptions, $L^{k+1}$ must be
rank $r$.

In order to show that $\tilde{L}^{k+1}$ is at most $2\mu$ incoherent, we apply
the following Lemma, which uses the Davis-Kahan inequality [davis-kahan].

######  Lemma 2.

If $X\in\mathbb{R}^{d_{1}\times d_{2}}$ (with $d_{1}\geq d_{2}$) is a rank
$r$, $\mu-$ incoherent matrix, and $\Delta\in\mathbb{R}^{d_{1}\times d_{2}}$
satisfies $||\Delta||_{2}\leq\frac{1}{2}\frac{\mu r}{d_{1}}\sigma_{r}(X)$,
then the top $r$ singular vectors of the matrix $X+\Delta$ form a $2\mu-$
incoherent tangent space.

Next, we bound $\Delta_{s}^{k}$ and showing $s^{k+1}\in\Omega$ in a similar
manner. We present the following Lemma, which mirrors Lemma 1.

######  Lemma 3.

Let $s^{*}\in\mathbb{R}^{d_{s}}$ be an sparse vector with support $\Omega$,
and let $\bar{S}$ be defined as

$\bar{s}=\underset{s\in\mathbb{R}^{d_{s}}}{\text{argmin
}}\;\lambda\phi(s)+||s-s^{*}+\delta||^{2}$ (11)

Where $\phi$ is an at most $\frac{1}{\lambda}$ weakly convex amenable
regualarizer, and $\delta\in\mathbb{R}^{d_{s}}$ satisfies
$||\delta||_{\infty}\leq\lambda$. Then, $\text{supp}(\bar{s})\subseteq\Omega$.
Furthermore, we can bound the difference $\Delta_{s}=\bar{s}-s^{*}$ as

$||\Delta_{s}||^{2}_{F}\leq||\mathcal{P}_{\Omega}\big{(}\delta\big{)}||_{F}+2\lambda\sqrt{|\Omega|}\phi^{\prime}\big{(}s_{min}-\lambda\big{)}.$
(12)

By the update equation for $\tilde{s}^{k+1}$, we have

$\displaystyle s^{*}-\tilde{s}^{k+1}=$ $\displaystyle
s^{*}-\big{(}s^{k}-\mathcal{A}_{s}^{*}\big{(}\mathcal{A}_{s}(s^{k})+\mathcal{A}_{L}(L^{k})-b)\big{)}\big{)}$
$\displaystyle=$ $\displaystyle
s^{*}-s^{k}+\frac{\tau_{s}}{d_{s}}\mathcal{A}_{s}^{*}\mathcal{A}_{s}s^{k}+\frac{\tau_{s}}{d_{s}}\mathcal{A}_{L}^{*}\mathcal{A}_{L}L^{k}-\frac{\tau_{s}}{d_{s}}A_{s}^{*}\big{(}\mathcal{A}_{s}(s^{*})+\mathcal{A}_{L}(L^{*})+\mathcal{E}\big{)}$
$\displaystyle=$
$\displaystyle\big{(}\Delta_{s}^{k}-\frac{\tau_{s}}{d_{s}}\mathcal{A}_{s}^{*}\mathcal{A}_{s}\Delta_{s}^{k}\big{)}-\frac{\tau_{s}}{d_{s}}\mathcal{A}_{s}^{*}\mathcal{A}_{L}\Delta_{L}^{k+1}-\frac{\tau_{s}}{d_{s}}\mathcal{A}_{s}^{*}\mathcal{E}$

By Lemma 3,

$\displaystyle||\Delta_{s}^{k+1}||_{F}^{2}\leq$
$\displaystyle||\mathcal{P}_{\Omega}\big{(}\Delta_{s}^{k}-\mathcal{A}_{s}^{*}\mathcal{A}_{s}\Delta_{s}^{k}\big{)}||_{F}^{2}+||\mathcal{P}_{\Omega}\mathcal{A}_{s}^{*}\mathcal{A}_{L}\Delta_{L}^{k+1}||_{F}^{2}$
$\displaystyle+||\mathcal{P}_{\Omega}\mathcal{A}_{s}^{*}\mathcal{E}||_{F}^{2}+\lambda_{s}|\Omega|\phi_{s}^{\prime}(s_{min}-\lambda_{s})$

Applying the RIP to the first term gives:

$||\mathcal{P}_{\Omega}\big{(}\Delta_{s}^{k}-\frac{\tau_{s}}{d_{s}}\mathcal{A}_{s}^{*}\mathcal{A}_{s}\Delta_{s}^{k}\big{)}||_{F}^{2}\leq\kappa_{s}||\Delta^{k}_{s}||_{F}^{2}$

And, applying the ROP to the second term gives us:

$||\frac{\tau_{s}}{d_{s}}\mathcal{P}_{\Omega}\mathcal{A}_{s}^{*}\mathcal{A}_{L}\Delta_{L}^{k+1}||_{F}^{2}=||\frac{\tau_{s}}{d_{s}}\mathcal{P}_{\Omega}\mathcal{A}_{s}^{*}\mathcal{A}_{L}\mathcal{P}_{T^{k+1}}\Delta_{L}^{k+1}||_{F}^{2}\leq\frac{\tau_{s}}{d_{s}}\kappa||\Delta_{L}^{k+1}||_{F}^{2}$

Combining these terms gives us the desired result.

###  3.4 Proofs of Supporting Lemmas

###### Proof of Lemma 3.

First, we will show that if $s^{*}_{i}=0$, then $\bar{s}_{i}=0$. For $i$ not
in $\Omega$,

$\bar{s}_{i}=\text{argmin}_{s}\;\;\lambda\phi(s)+\frac{1}{2}(s-\delta_{i})^{2}$

The sub-gradient of the objective function evaluated at zero is
$\\{\delta_{i}+g\;\big{|}\;|g|\leq\lambda\\}$. By our assumption that
$\lambda\geq\delta_{i}$, 0 must be in this set, and so $\bar{s}_{i}=0$ is a
stationary point. Because the objective function is strongly convex by
assumption, this is a global minimizer.

Next, consider if $i\in\Omega$. By first order necessary conditions for
optimality,

$\displaystyle\bar{s}_{i}=s_{i}^{*}+\delta_{i}-\phi^{\prime}(\bar{s}_{i})\leq
s_{i}^{*}-\delta_{i}-\phi^{\prime}(s_{ij}^{*}+\delta_{i})$

This gives the bound:

$|\bar{s}_{i}-s_{i}^{*}|\leq|\delta_{i}|+\phi^{\prime}(s_{min}-\lambda)$

Combining these facts,

$\displaystyle||\bar{s}-s^{*}||^{2}=||\mathcal{P}_{\Omega}(\bar{s}-s^{*})||^{2}\leq||\mathcal{P}_{\Omega}(\delta)||^{2}+|\Omega|\phi^{\prime}(s_{min}-\lambda)^{2}$

∎

###### Proof of Lemma 1.

Let

$\displaystyle f(L)=\lambda\Phi(L)+\frac{1}{2}||L-L^{*}+\delta||_{F}^{2}$

be the function that $\bar{L}$ is the global minimizer of, that is,
$0\in\partial f(\bar{L})$. The subdifferential of $f$ at $L^{*}$ is as
follows:

$\big{\\{}\lambda\big{(}U\Sigma^{\phi}V^{T}+W\big{)}+\delta\>\bigg{|}\>W\in
T^{\perp},||W||_{2}\leq 1\big{\\}}$

where $L^{*}=U\Sigma V^{T}$ and

$\Sigma^{\phi}=\text{diag}(\phi^{\prime}(\sigma_{1}(L^{*})),\phi^{\prime}(\sigma_{2}(L^{*})),\ldots,\phi^{\prime}(\sigma_{r}(L^{*}))).$

Because $\Phi$ is $\frac{1}{2\lambda}$-weakly convex, $f$ is
$\frac{1}{2}$-strongly convex, which gives us:

$\displaystyle||L^{*}-\bar{L}||_{F}\leq$ $\displaystyle\inf 2||\partial
f(L^{*})||_{F}.$

Consider the subgradient given by
$W=-\frac{1}{\lambda}\mathcal{P}_{\mathcal{T}^{\perp}}(\delta)$. Note the that
$||W||_{2}\leq 1$ as we assume $\lambda\geq||\delta||_{2}$.

$\displaystyle||L^{*}-\bar{L}||_{F}\leq$ $\displaystyle 2||\lambda
U\Sigma^{\phi}V^{T}+\delta-\mathcal{P}_{\mathcal{T}^{\perp}}(\delta)||_{F}$
$\displaystyle\leq$ $\displaystyle
2\lambda||U\Sigma^{\phi}V^{T}||_{F}+2||\delta-\mathcal{P}_{\mathcal{T}^{\perp}}(\delta)||_{F}$
$\displaystyle\leq$ $\displaystyle
2\lambda\sqrt{r}\phi^{\prime}(\sigma_{r}(L^{*}))+2||\mathcal{P}_{\mathcal{T}}(\delta)||_{F}$

The second inequality is the triangle inequality, and the third uses the fact
that $||\Sigma^{\phi}||_{2}=\phi^{\prime}(\sigma_{r}(L^{*}))$ ∎

###### Proof of Lemma 2.

Let ${U}\in\mathbb{R}^{d_{1}\times r}$ and
$\tilde{U}\in\mathbb{R}^{d_{1}\times r}$ denote the (top $r$) left singular
vectors of $X$ and $X+\Delta$ respectively. By the Davis Kahan theorem,

$\text{dist}\big{(}U,\tilde{U}\big{)}\leq\frac{||\Delta||_{2}}{\sigma_{r}(X)+||\Delta||_{2}}\leq\frac{1}{2}\frac{\mu
r}{d_{1}}$

where the second inequality uses the fact that $||\Delta||_{2}\leq\frac{\mu
r}{2d_{1}}\sigma_{r}(X)$ Let $u_{i}$ and $\tilde{u}_{i}$ be the $i^{th}$ left
singular vectors of $X$ and $X+\Delta$ respectively, and let $\theta_{i}$ be
the angle between the vectors.

$\max(\tilde{u}_{i})\leq\max(u_{i})+2\sin(\theta_{i})\leq\frac{\mu
r}{d_{1}}+\frac{\mu r}{d_{1}}=\frac{2\mu r}{d_{1}}$

where we use the fact that

$\text{dist}\big{(}U,\tilde{U}\big{)}=\text{max}_{i}\>\sin(\theta_{i})$

So, the rank $r$ approximation of $X+\Delta$ is $2\mu$ incoherent. ∎

##  4 Results for Specific Models

In this section, we use Theorem 1 to analyze an application of the alternating
proximal gradient descent algorithm to the problems of matrix completion and
RPCA.

###  4.1 Matrix Completion

We start by considering the problem of matrix completion. Here, we have a
sparse set of observed entries of $L$, $\mathcal{A}_{L}=\mathcal{A}_{\Omega}$,
and we do not consider a sparse vector $s$ (i.e. $A_{s}=0$).

We present a version of the RIP for the sampling operator from [29].

######  Lemma 4 ([29]).

Let $\Omega$ be a set of $n$ entries of
$\\{1,\ldots,d_{1}\\}\times\\{1,\ldots,d_{2}\\}$ drawn independently at random
with uniform probability, with $n>64\mu r(d_{1}+d_{2})\log(d_{2})$. Then, with
probability at least $1-2d_{2}^{-2}$,

$\frac{5}{6}||X||_{F}^{2}\leq\frac{d_{1}d_{2}}{n}||\mathcal{A}_{\Omega}(X)||^{2}\leq\frac{7}{6}||X||_{F}^{2}$

for any rank $r$, $\mu-$incoherent matrix $X$.

Additionally, we will assume that the additive noise $\mathcal{E}$ has entries
that are mean zero i.i.d. variables. The effect the noise has on the estimator
is reduced due to the fact that very little of
$\mathcal{A}^{*}_{\Omega}(\mathcal{E})$ will lie in the tangent space
$\mathcal{T}$. To formalize this intuition, we cite the following lemma from
[38].

######  Lemma 5.

Assume that the entries of $\Omega$ are chosen uniformly at random from
$\\{1,\ldots,d_{1}\\}\times\\{1,\ldots,d_{2}\\}$, with
$d=\text{max}(d_{1},d_{2})$, and the entries of $\mathcal{E}$ are mean zero
i.i.d. random variables with variance $\nu^{2}$. For some universal constants
$C_{1}$ and $C_{2}$, the following hold:

$||\mathcal{A}^{*}_{\Omega}(\mathcal{E})||_{2}\leq
C_{1}\nu\sqrt{pdlog(d)}\;\;\text{ and
}\;\;||\mathcal{A}^{*}_{\Omega}(\mathcal{E})||_{\infty}\leq
C_{2}\nu\sqrt{plog(d)}$

We now present an error bound of for a stationary point of Algorithm 1 applied
to matrix completion.

######  Theorem 2.

Let $b=\mathcal{A}_{\Omega}X^{*}+\mathcal{E}$ for a rank $r$, $\mu-$
incoherent matrix $X^{*}\in\mathbb{R}^{d_{1}\times d_{2}}$, and the entries of
$\mathcal{E}$ are mean zero i.i.d. random variables with variance $\nu^{2}$.
There exists universal constants $C_{1}$ and $C_{2}$ such that under the same
assumptions as Lemma 4, if $\lambda>C_{1}\nu\sqrt{pd\text{ log}(d)}$, then the
iterates of Algorithm 1 linearly converge to a point $\bar{X}$ such that
$||X^{*}-\bar{X}||_{F}^{2}$ is less than:

$C_{2}\frac{d_{1}d_{2}}{n}\big{(}\underbrace{r\nu^{2}d\text{
log}(d)}_{\begin{array}[]{c}\text{ optimal}\\\ \text{error
rate}\end{array}}+\underbrace{\lambda
r\phi^{\prime}(\sigma_{r}(X^{*}))}_{\text{bias term}}\big{)}$

with convergence rate $\frac{1}{6}\frac{d_{1}d_{2}}{n}$.

The two terms of the error bound account for the optimal error rate and a bias
term. The optimal error rate is the error bound if we know the tangent space
of $X^{*}$ a priori, that is, the difference between $X^{*}$ and the solution
to the optimization problem

$\underset{X\in\mathcal{T}}{\text{min
}}\;||\mathcal{A}_{\Omega}(X)-b||_{F}^{2}.$

The oracle rate is further discussed in [7] and [26].

###### Proof.

We can apply Theorem 1 with $s^{k}=s^{*}=0$, $\tau_{L}=1$, and
$\kappa_{L}=\frac{1}{6}$ (by Lemma 3) to get the bound

$||\Delta^{k+1}||_{F}^{2}\leq\frac{1}{6}||\Delta^{k}||_{F}^{2}+\big{(}\frac{d_{1}d_{2}}{n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{\Omega}^{*}\mathcal{E}||_{F}^{2}+\lambda
r\phi^{\prime}(\sigma_{r}(X^{*})-2\lambda)\big{)}$

Initializing with $X^{0}=0$, we have the error at each iteration as follows:

$||\Delta^{k}||^{2}\leq\big{(}\frac{1}{6}\big{)}^{k}||X^{*}||_{F}^{2}+\frac{6}{5}\big{(}1-\big{(}\frac{1}{6}\big{)}^{k}\big{)}\big{(}\frac{d_{1}d_{2}}{n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{\Omega}^{*}\mathcal{E}||_{F}^{2}+\lambda
r\phi^{\prime}(\sigma_{r}(X^{*})-2\lambda)\big{)}$

Taking the limit as $k\rightarrow\infty$, and applying the bounds from Lemma 5
gives the desired result. ∎

Perhaps counter-intuitively, a choice of step size, $\tau$, that minimizes the
loss function (i.e. the steepest descent step size) is not always the step
size that leads to the fastest convergence rate. To see this, we compare the
error bound at iteration $k$ in both cases. The steepest descent step size
($\tau=n$) would give

$\displaystyle L^{*}-\tilde{L}^{k}$
$\displaystyle=L^{*}-L^{k-1}-\mathcal{P}_{\Omega}(L^{*}-L^{k-1})=P_{\Omega^{C}}(\Delta_{L})$

where $\Omega^{C}$ is the set of indices not in $\Omega$. The stepsize
informed by the RIP ($\tau=\frac{d_{1}d_{2}}{n}$) gives

$\displaystyle L^{*}-\tilde{L}^{k}$
$\displaystyle=L^{*}-L^{k-1}-\frac{d_{1}d_{2}}{n}\mathcal{P}_{\Omega}(L^{*}-L^{k-1})=\Delta_{L}-\frac{d_{1}d_{2}}{n}P_{\Omega}(\Delta_{L}).$

While the $\tau=1$ gives a significantly smaller error when simply comparing
$L^{*}-\tilde{L}^{k}$, the error bound for $L^{*}-L^{k}$ comes from projecting
$L^{*}-\tilde{L}^{k}$ onto $\mathcal{T}$. Without further information, the
best bound we can get when using $\tau=n$ would be

$||L^{*}-L^{k}||_{F}^{2}\leq
1-\frac{(1-\kappa)n}{d_{1}d_{2}}||L^{*}-L^{k-1}||_{F}^{2}.$

This convergence rate approaches 1 asymptotically when we consider the
information theoretic minimum number of measurements for large $d_{1}$ and
$d_{2}$. However, when $\tau=\frac{d_{1}d_{2}}{n}$, the convergence rate
remains constant:

$\displaystyle||L^{*}-L^{k}||_{F}^{2}$
$\displaystyle\leq||\Delta_{L}-\frac{d_{1}d_{2}}{n}\mathcal{P}_{\mathcal{T}}\mathcal{P}_{\Omega}\mathcal{P}_{\mathcal{T}}(\Delta_{L})||_{F}^{2}\leq\kappa||\Delta_{L}||_{F}^{2}$

The first inequality uses the fact that $\Delta_{L}\in\mathcal{T}$ and the
second uses the RIP.

###  4.2 Robust PCA

Next, we will use Theorem 1 to analyze the APGD algorithm applied to the
problem of RPCA. Specifically, we are interested in the special case of (1)
where the $A_{s}=I_{n}$, and $\mathcal{A}_{L}=\mathcal{A}_{\Omega^{obs}}$.

In order for RPCA to be possible, we need the nonzero entries of
$\mathcal{A}_{\Omega^{obs}}^{*}s^{*}$ to be sufficiently well-distributed
throughout the rows and columns – if the sparse corruptions affected the same
row or column of $L^{*}$, then this would also be a low-rank perturbation and
thus be impossible to separate from $L^{*}$ without further information. So,
we will assume that $\mathcal{A}_{L}^{*}s^{*}$ is $\alpha-$sparse, defined as
follows.

######  Definition 5.

The matrix $S$ is $\alpha$-sparse for $0<\alpha<1$ if the proportion of
nonzero entries in any row or column is less than $\alpha$. That is,

$||S_{i:}||_{0}\leq\alpha d_{1}\>,||S_{:j}||_{0}\leq\alpha d_{2}\;\forall i,j$
(13)

In order to verify that the ROP property holds, we present the following
Lemma:

######  Lemma 6.

Let $\mathcal{T}$ be a rank $r$, $\mu-$incoherent tangent space, and let
$\Omega$ be an $\alpha-$sparse subspace. Then,

$||\mathcal{P}_{\mathcal{T}}\mathcal{P}_{\Omega}||\leq 2\alpha\mu
r,\;||\mathcal{P}_{\Omega}\mathcal{P}_{\mathcal{T}}||\leq 2\alpha\mu r$ (14)

###### Proof.

By the triangle inequality, for a matrix $S\in\Omega$,

$\displaystyle||\mathcal{P}_{\mathcal{T}}(S)||_{F}^{2}\leq||UU^{T}S||^{2}_{F}+||SVV^{T}||_{F}^{2}+||UU^{T}SVV^{T}||_{F}^{2}$

Because $U$ is an orthonormal matrix,

$\displaystyle||UU^{T}S||^{2}_{F}=$
$\displaystyle\text{trace}(UU^{T}SS^{T}UU^{T})=\text{trace}(U^{T}UU^{T}SS^{T}U)$
$\displaystyle=$ $\displaystyle\text{trace}(U^{T}SS^{T}U)=||U^{T}S||_{F}^{2}$

We now expand this norm and use the incoherence property to obtain the desired
result:

$\displaystyle||U^{T}S||_{F}^{2}=$
$\displaystyle\sum_{k=1}^{r}\sum_{i=1}^{d_{1}}\langle
U_{k},S_{i}\rangle^{2}\leq\sum_{k=1}^{r}\sum_{i=1}^{d_{1}}\big{(}\sum_{j\in\Omega_{i}}U_{kj}^{2}\big{)}||S_{i}||^{2}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{d_{1}}||S_{i}||^{2}\sum_{j\in\Omega_{i}}||U_{j}||^{2}\leq\sum_{i=1}^{d_{1}}||S_{i}||^{2}\alpha\mu
r\leq\alpha\mu r||S||_{F}^{2}$

∎

We can now give a bound for the stationary point of the APGD algorithm applied
to RPCA.

######  Theorem 3.

Let $b=\mathcal{A}_{\Omega^{obs}}(L^{*})+s^{*}+\mathcal{E}$ for a rank $r$,
$\mu-$ incoherent matrix $L^{*}\in\mathbb{R}^{d_{1}\times d_{2}}$, a sparse
vector $s^{*}\in\mathbb{R}^{n}$ with $||s^{*}||_{\infty}\leq
2||L^{*}||_{\infty}$, and a vector $\mathcal{E}\in\mathbb{R}^{n}$ whose
entries are independent Gaussian variables with mean 0 and standard deviation
$\sigma$. Under the same assumptions on $\Omega^{obs}$ as Lemma 4, and
assuming that $\mathcal{A}_{\Omega}^{*}s^{*}$ is $\alpha-$ sparse and
$\alpha\mu r\leq\frac{1}{64}$, if
$\lambda_{L}\geq\frac{1}{6}+||\mathcal{A}_{\Omega}^{*}\mathcal{E}||_{2}$ and
$\lambda_{s}\geq\frac{\mu r}{d_{1}d_{2}}+||\mathcal{E}||_{\infty}$, then the
iterates of Algorithm 1 linearly converge to a point $\bar{L},\bar{s}$
satisfying

$\displaystyle||\Delta_{L}||_{F}^{2}\leq
C_{1}\frac{d_{1}d_{2}}{n}r\nu^{2}d\text{
log}(d),\;\;\>||\Delta_{S}||_{F}^{2}\leq
C_{2}\frac{d_{1}d_{2}}{n}r\nu^{2}d\text{ log}(d)$

with convergence rate $\frac{1}{6}$.

###### Proof.

In order to apply Lemma 4, we first show the following properties about the
parameters $\lambda_{L}$ and $\lambda_{s}$.

$\displaystyle\lambda_{L}$
$\displaystyle\geq||\mathcal{A}_{L}^{*}\mathcal{E}||_{2}+||\mathcal{A}_{L}^{*}s^{*}||_{2}$
$\displaystyle\geq||\mathcal{A}_{L}^{*}\mathcal{E}||_{2}+||\mathcal{A}_{L}^{*}\delta_{s}^{k}||_{2}\;\forall
k$

where the first inequality comes from the assumption on $\lambda_{L}$ and the
second from the fact that the two norm of $s^{k}$ is decreasing.

$\displaystyle\lambda_{s}\geq$
$\displaystyle||\mathcal{E}||_{\infty}+\frac{\mu r}{d_{1}d_{2}}$
$\displaystyle\geq$ $\displaystyle||\mathcal{E}||_{\infty}+||L^{*}||_{\infty}$
$\displaystyle\geq$
$\displaystyle||\mathcal{E}||_{\infty}+||\Delta_{L}^{k}||_{\infty}\;\forall k$

where the first inequality is our assumption on $\lambda_{s}$, the second
follows from incoherence of $L^{*}$, and the last from the fact that

$||\Delta_{L}^{k}||_{\infty}\leq||\Delta_{L}^{0}||_{\infty}=||L^{*}||_{\infty}.$

In order to show the bias term is 0, we use our assumption that
$\lambda_{L}\leq\frac{1}{4}\sigma_{r}(L^{*})$ and the fact that the MCP
regularizer can satisfy $\phi_{L}^{\prime}(3\lambda_{L})=0$ while still
maintaining $\frac{1}{2\lambda_{L}}$-weak convexity. Likewise, when we assume
$\lambda_{s}\leq\frac{1}{4}s_{min}$, we can choose $\phi_{s}$ such that
$\phi^{\prime}_{s}(s_{min}-\lambda_{s})\leq\phi^{\prime}_{s}(3\lambda_{s})=0$
while still maintaining $\frac{1}{2\lambda_{s}}$ weak convexity.

We can now apply the first part of Theorem 1 to obtain the bound:

$\displaystyle||\Delta^{k+1}_{L}||_{F}^{2}$
$\displaystyle\leq\frac{1}{6}||\Delta^{k}_{L}||_{F}^{2}+\frac{1}{64}||\Delta^{k}_{s}||_{F}^{2}+\frac{d_{1}d_{2}}{n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}\mathcal{E}||_{F}^{2}$
$\displaystyle||\Delta_{S}^{k+1}||_{F}^{2}$
$\displaystyle\leq\frac{1}{64}||\Delta_{L}^{k+1}||_{F}^{2}+||\mathcal{P}_{\Omega}\mathcal{E}||_{F}^{2}$
$\displaystyle\leq\frac{1}{384}||\Delta_{L}^{k}||_{F}^{2}+\frac{1}{4096}||\Delta_{S}^{k}||_{F}^{2}+\frac{d_{1}d_{2}}{64n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}\mathcal{E}||_{F}^{2}+||\mathcal{P}_{\Omega}\mathcal{E}||_{F}^{2}$

The limit point of this sequence gives:

$\displaystyle||\Delta_{L}||_{F}^{2}$
$\displaystyle\leq\frac{6}{5}\frac{d_{1}d_{2}}{n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}\mathcal{E}||_{F}^{2}$
$\displaystyle||\Delta_{S}||_{F}^{2}$
$\displaystyle\leq\frac{d_{1}d_{2}}{32n}||\mathcal{P}_{\mathcal{T}}\mathcal{A}_{L}^{*}\mathcal{E}||_{F}^{2}+\frac{6}{5}||\mathcal{P}_{\Omega}\mathcal{E}||_{F}^{2}$

From here, we can apply the bounds from Lemma 5 to get the desired results. ∎

The bound presented in Theorem 3 are a major improvement upon previous bounds,
such as the ones presented in Agarwal, Negahban, and Wainwright [2]. For the
estimator (2), they require that

$\lambda_{s}\geq 4(||L^{*}||_{\infty}+||E||_{\infty}).$

However, because the sparse matrix $S$ obtained from their method will be zero
whenever the true value is less than $\lambda_{s}$, their result has an
implicit assumption that the nonzero values of $S^{*}$ are four times larger
than any value in $L^{*}+E$. If this were true, then identifying entries of
the measured matrix $M$ could be accomplished by identifying the entries with
the largest absolute value. Our result imposed the much less strict assumption
that the nonzero entries of $S^{*}$ are larger than $||E||_{\infty}$ plus a
small constant (an assumption that is necessary for $S^{*}$ and $E$ to be
separable).

##  5 Numerical Results

We implemented Algorithm 1 in Matlab R2020a, for which the code is available
at GitHub. All results in this section are obtained with the Matlab version in
order to accurately compare to other algorithms which are only available in
Matlab, and are run on a Windows 10 desktop with an AMD Phenom 3.40 GHz
processor and 8 Gb of RAM.

###  5.1 Matrix Completion

Table 2: Comparison of four different matrix completion algorithms on randomly generated low rank matrices and common recommendation data sets. The algorithm LMaFit reconstructs a matrix of a given rank $k$. The table shows the results when the algorithm is given the exact rank($k=r$) and an incorrect rank ($k=2r$). Time is given in seconds  $d_{1}$ | 1000 | 1000 | 5000 | 5000
---|---|---|---|---
$d_{2}$ | 500 | 500 | 1000 | 1000
$r$ | 5 | 5 | 10 | 10
$n/d_{1}d_{2}$ | 0.3 | 0.1 | 0.2 | 0.05
std($E$) | 0.1 | 0.02 | 0.1 | 0.02
| RFNE | T | RFNE | T | RFNE | T | RFNE | T
APGD | 3.28e-4 | 0.7 | 2.90e-4 | 1.1 | 1.69e-4 | 10 | 1.96e-4 | 29
FaNCL | 3.28e-4 | 1.8 | 2.90e-4 | 4.3 | 1.69e-4 | 31 | 2.49E-4 | 57
IALM | 3.28e-4 | 2.6 | 2.92e-4 | 2.8 | 1.72e-4 | 32 | 1.99e-4 | 27
LMaFit (2$r$) | 4.68e2 | 21 | 1.20e3 | 6.2 | 4.99e2 | 197 | 1.60e3 | 40
LMaFit($r$) | 3.28e-4 | 0.4 | 2.90e-4 | 0.3 | 1.69e-4 | 3.9 | 1.96e-4 | 3.9
Table 3: Comparison of four different matrix completion algorithms on two common recommender system benchmarking datasets.  | Jester Dataset | ML 1M
---|---|---
Number of Users ($d_{1}$) | 24983 | 6040
Number of Items ($d_{2}$) | 100 | 3952
Percentage of Entries Observed ($n/d_{1}d_{2}$) | 0.58 | 0.034
| NMAE | Time (s) | NMAE | Time (s)
APGD | 0.159 | 21 | 0.172 | 172
FaNCL | 0.183 | 42 | 0.200 | 42
IALM | 0.163 | 77 | 0.183 | 216
LMaFit | 0.168 | 7.4 | 9.174 | 44

We compare our method for matrix completion to another method utilizing
nonconvex regularizer from [36] (FaNCL), along with a method to minimize the
nuclear norm IALM, from [19], and a rank constrained method, LMaFit [35]. The
results are shown in Table 2.

We start by comparing the performance of the methods on randomly generated low
rank matrices of varying size, rank, percentage of observed entries, and
standard deviation of the noise in the measurements (shown relative to the
mean absolute value of the low rank matrix). In each of the cases, our method
performs exactly as well as FaNCL and LMaFit when the correct rank is given.
IALM performs equally well in the first case, and slightly worse than the
remaining three cases due to the fact that the nuclear norm biases the result
towards zero.

Next, we show the results on common data sets for recommendation systems, the
Jester data set [14] and MovieLens 1M [1]. For each of these two data sets, we
partition the observations into five folds, fit a low rank model to four of
the folds and calculate the accuracy on the remaining fold. We repeat this for
each of the five folds and present the average normalized mean absolute error.
In both cases, our method outperforms the other three algorithms we compare
to.

###  5.2 Robust PCA

We compare our method to several other prominent RPCA methods, including
LMaFit [31], AltProj [27], RPCA-GD [37], and IALM [19]. We compared with many
other methods included in the LRSLibrary [32], however we only include results
from the aforementioned algorithms as they gave the most accurate results for
matrices with a significant amount with noise, a test case we emphasise in
this section.

It is worth noting that, while our algorithm does not require an estimate of
the rank of $L^{*}$ or the sparsity of $S^{*}$ a priori, RPCA-GD requires
both, and LMaFit and AltProj require an estimate of the rank. However, we
found that LMaFit and AltProj still perform very well when this estimate is
unreliable, as LMaFit includes a rank-estimation scheme and AltProj starts by
performing a rank 1 projection, and increases the rank up until the estimate
given. We provide all methods with an upper bound on the rank equal to twice
the rank of $L^{*}$ and an upper bound on the number of corrupted entries
equal to twice that of $S^{*}$.

Table 4: Accuracy and run-time for RPCA on four videos of fish. | Marine Snow | Aquaculture | Caustics | Two Fish | Fish Swarm
---|---|---|---|---|---
| Acc | Time | Acc | Time | Acc | Time | Acc | Time | Acc | Time
APGD | 0.92 | 87 | 0.92 | 592 | 0.88 | 459 | 0.90 | 467 | 0.75 | 473
LMaFit | 0.59 | 258 | 0.65 | 316 | 0.7 | 273 | 0.62 | 291 | 0.55 | 310
AltProj | 0.84 | 143 | 0.62 | 109 | 0.86 | 83 | 0.76 | 109 | 0.62 | 130
RPCA-GD | 0.92 | 94 | 0.72 | 77 | 0.87 | 76 | 0.96 | 96 | 0.66 | 93
IALM | 0.57 | 379 | 0.56 | 286 | 0.56 | 288 | 0.56 | 342 | 0.56 | 305
Figure 1: Frames from the Marine Snow video. The first row of the first column
shows the original frame, with the image segmentation for that frame below it.
The remaining three columns show the background and foreground obtained by
three different methods.  Figure 2: Frames from the Small Aquaculture video.
Figure 3: Frames from the Small Aquaculture video. Figure 4: Frames from the
Small Aquaculture video. Figure 5: Frames from the Small Aquaculture video.

The most commonly used test case in the literature for RPCA is the task of
separating the background and foreground of a video. In this scenario, each
frame of the video is represented as a column vector in $M$. If the background
is static (or, at least, in some way repetitive), then we can expect the
background of each frame to be represented as a low rank matrix, and the
foreground as a sparse matrix. See [6] for further details.

We evaluated our algorithm on the Underwater Change Detection dataset [28].
Out of the 1100 frames in each video, the ground truth image segmentation is
included for last 100 frames. This allows use to present an objective and
accurate metric of how well each method is able to identify the foreground of
the image.

In Table 4, we present the runtime and the accuracy of determining which
pixels contain a fish for our method compared to the four previously mentioned
approaches. To calculate the accuracy, we average of the true positive ratio
and true negative ratio. In three of the five videos, our method achieves the
highest accuracy, whereas in the other two the RPCA-GD algorithm preforms
slightly better. The recovered background and foreground for our method, RPCA-
GD and AltProj are shown in Figures 1, 2, 3, 4, 5, along with the original
frame and segmented image.

##  6 Conclusions

We have shown a novel convergence analysis of the alternating proximal
gradient descent algorithm applied to the problems of matrix completion and
RPCA with nonconvex regularizers, and bound the difference from the ground
truth low rank matrix and sparse vector. Future work on the topic could
include extending our analysis to data that lies on more complicated,
nonlinear manifolds.

## References

  * [1] Movielens.  https://grouplens.org/datasets/movielens/.  Accessed: 2019-11-21.
  * [2] A. Agarwal, S. Negahban, and M. J. Wainwright, Noisy matrix decomposition via convex relaxation: Optimal rates in high dimensions, Ann. Statist., 40 (2012), pp. 1171–1197, https://doi.org/10.1214/12-AOS1000, https://doi.org/10.1214/12-AOS1000.
  * [3] J.-F. Cai, E. J. Candes, and Z. Shen, A singular value thresholding algorithm for matrix completion, 2008, https://arxiv.org/abs/0810.3286.
  * [4] E. Candès and B. Recht, Exact matrix completion via convex optimization, Commun. ACM, 55 (2012), p. 111–119, https://doi.org/10.1145/2184319.2184343, https://doi.org/10.1145/2184319.2184343.
  * [5] E. Candes and T. Tao, Decoding by linear programming, IEEE Transactions on Information Theory, 51 (2005), pp. 4203–4215, https://doi.org/10.1109/TIT.2005.858979.
  * [6] E. J. Candès, X. Li, Y. Ma, and J. Wright, Robust principal component analysis?, J. ACM, 58 (2011), https://doi.org/10.1145/1970392.1970395, https://doi.org/10.1145/1970392.1970395.
  * [7] E. J. Candes and Y. Plan, Matrix completion with noise, Proceedings of the IEEE, 98 (2010), pp. 925–936, https://doi.org/10.1109/JPROC.2009.2035722.
  * [8] E. J. Candes and T. Tao, The power of convex relaxation: Near-optimal matrix completion, IEEE Transactions on Information Theory, 56 (2010), pp. 2053–2080, https://doi.org/10.1109/TIT.2010.2044061.
  * [9] E. J. Candes, M. B. Wakin, and S. P. Boyd, Enhancing sparsity by reweighted $l$ 1 minimization, Journal of Fourier analysis and applications, 14 (2008), pp. 877–905.
  * [10] V. Chandrasekaran, S. Sanghavi, P. A. Parrilo, and A. S. Willsky, Rank-sparsity incoherence for matrix decomposition, SIAM Journal on Optimization, 21 (2011), pp. 572–596, https://doi.org/10.1137/090761793, https://doi.org/10.1137/090761793, https://arxiv.org/abs/https://doi.org/10.1137/090761793.
  * [11] R. Chartrand, Nonconvex splitting for regularized low-rank + sparse decomposition, IEEE Transactions on Signal Processing, 60 (2012), pp. 5810–5819, https://doi.org/10.1109/TSP.2012.2208955.
  * [12] R. Ge, C. Jin, and Y. Zheng, No spurious local minima in nonconvex low rank problems: A unified geometric analysis, vol. 70 of Proceedings of Machine Learning Research, International Convention Centre, Sydney, Australia, 06–11 Aug 2017, PMLR, pp. 1233–1242, http://proceedings.mlr.press/v70/ge17a.html.
  * [13] R. Ge, J. D. Lee, and T. Ma, Matrix completion has no spurious local minimum, in Advances in Neural Information Processing Systems, D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, eds., vol. 29, Curran Associates, Inc., 2016, pp. 2973–2981, https://proceedings.neurips.cc/paper/2016/file/7fb8ceb3bd59c7956b1df66729296a4c-Paper.pdf.
  * [14] K. Goldberg, T. Roeder, D. Gupta, and C. Perkins, Eigentaste: A constant time collaborative filtering algorithm, Inf. Retr., 4 (2001), p. 133–151, https://doi.org/10.1023/A:1011419012209, https://doi.org/10.1023/A:1011419012209.
  * [15] T. A. Horn and C. R. Johnson, Topics in Matrix Analysis, Cambridge, 1991\.
  * [16] P. Jain, R. Meka, and I. Dhillon, Guaranteed rank minimization via singular value projection, in Advances in Neural Information Processing Systems, J. Lafferty, C. Williams, J. Shawe-Taylor, R. Zemel, and A. Culotta, eds., vol. 23, Curran Associates, Inc., 2010, pp. 937–945, https://proceedings.neurips.cc/paper/2010/file/08d98638c6fcd194a4b1e6992063e944-Paper.pdf.
  * [17] P. Jain, P. Netrapalli, and S. Sanghavi, Low-rank matrix completion using alternating minimization, in Proceedings of the Forty-Fifth Annual ACM Symposium on Theory of Computing, STOC ’13, New York, NY, USA, 2013, Association for Computing Machinery, p. 665–674, https://doi.org/10.1145/2488608.2488693, https://doi.org/10.1145/2488608.2488693.
  * [18] Z. Kang, C. Peng, and Q. Cheng, Robust pca via nonconvex rank approximation, in 2015 IEEE International Conference on Data Mining, 2015, pp. 211–220, https://doi.org/10.1109/ICDM.2015.15.
  * [19] Z. Lin, M. Chen, and Y. Ma, The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices, arXiv preprint arXiv:1009.5055, (2010).
  * [20] P.-L. Loh and M. J. Wainwright, Regularized m-estimators with nonconvexity: Statistical and algorithmic theory for local optima, Journal of Machine Learning Research, 16 (2015), pp. 559–616, http://jmlr.org/papers/v16/loh15a.html.
  * [21] P.-L. Loh and M. J. Wainwright, Support recovery without incoherence: A case for nonconvex regularization, Ann. Statist., 45 (2017), pp. 2455–2482, https://doi.org/10.1214/16-AOS1530, https://doi.org/10.1214/16-AOS1530.
  * [22] C. Lu, J. Tang, S. Yan, and Z. Lin, Nonconvex nonsmooth low rank minimization via iteratively reweighted nuclear norm, IEEE Transactions on Image Processing, 25 (2016), pp. 829–839, https://doi.org/10.1109/TIP.2015.2511584.
  * [23] C. Lu, C. Zhu, C. Xu, S. Yan, and Z. Lin, Generalized singular value thresholding, Procedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence, (2015), pp. 1805–1811.
  * [24] K. Mohan and M. Fazel, Iterative reweighted algorithms for matrix rank minimization, J. Mach. Learn. Res., 13 (2012), p. 3441–3473.
  * [25] S. Negahban and M. J. Wainwright, Estimation of (near) low-rank matrices with noise and high-dimensional scaling, The Annals of Statistics, (2011), pp. 1069–1097.
  * [26] S. Negahban and M. J. Wainwright, Restricted strong convexity and weighted matrix completion: Optimal bounds with noise, The Journal of Machine Learning Research, 13 (2012), pp. 1665–1697.
  * [27] P. Netrapalli, U. Niranjan, S. Sanghavi, A. Anandkumar, and P. Jain, Non-convex robust pca, in Advances in Neural Information Processing Systems, 2014, pp. 1107–1115.
  * [28] M. Radolko, F. Farhadifard, and U. F. von Lukas, Dataset on underwater change detection, in OCEANS 2016 MTS/IEEE Monterey, 2016, pp. 1–8, https://doi.org/10.1109/OCEANS.2016.7761129.
  * [29] B. Recht, A simpler approach to matrix completion, J. Mach. Learn. Res., 12 (2011), p. 3413–3430.
  * [30] A. Sagan and J. E. Mitchell, Low-rank factorization for rank minimization with nonconvex regularizers, 2021.
  * [31] Y. Shen, Z. Wen, and Y. Zhang, Augmented lagrangian alternating direction method for matrix separation based on low-rank factorization, Optimization Methods and Software, 29 (2014), pp. 239–263.
  * [32] A. Sobral, T. Bouwmans, and E.-h. Zahzah, Lrslibrary: Low-rank and sparse tools for background modeling and subtraction in videos, in Robust Low-Rank and Sparse Matrix Decomposition: Applications in Image and Video Processing, CRC Press, Taylor and Francis Group., 2015.
  * [33] B. Vandereycken, Low-rank matrix completion by riemannian optimization, SIAM Journal on Optimization, 23 (2013), pp. 1214–1236.
  * [34] Z. Wang, H. Liu, and T. Zhang, Optimal computational and statistical rates of convergence for sparse nonconvex learning problems, Ann. Statist., 42 (2014), pp. 2164–2201, https://doi.org/10.1214/14-AOS1238, https://doi.org/10.1214/14-AOS1238.
  * [35] Z. Wen, W. Yin, and Y. Zhang, Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm, Mathematical Programming Computation, 4 (2012), pp. 333–361.
  * [36] Q. Yao, J. T.-Y. Kwok, and B. Han, Efficient nonconvex regularized tensor completion with structure-aware proximal iterations, vol. 97 of Proceedings of Machine Learning Research, Long Beach, California, USA, 09–15 Jun 2019, PMLR, pp. 7035–7044, http://proceedings.mlr.press/v97/yao19a.html.
  * [37] X. Yi, D. Park, Y. Chen, and C. Caramanis, Fast algorithms for robust pca via gradient descent, in Advances in Neural Information Processing Systems, D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, eds., vol. 29, Curran Associates, Inc., 2016, pp. 4152–4160, https://proceedings.neurips.cc/paper/2016/file/b5f1e8fb36cd7fbeb7988e8639ac79e9-Paper.pdf.
  * [38] X. Zhang, L. Wang, and Q. Gu, A unified framework for nonconvex low-rank plus sparse matrix recovery, vol. 84 of Proceedings of Machine Learning Research, Playa Blanca, Lanzarote, Canary Islands, 09–11 Apr 2018, PMLR, pp. 1097–1107, http://proceedings.mlr.press/v84/zhang18c.html.