text:
\section{Introduction}

\IEEEPARstart{P}{rincipal} component analysis (PCA) is a fundamental tool in statistics and data science. It obtains a low-dimensional expression for high-dimensional data in an $\ell_2$ sense. However, it is known that the classical PCA is sensitive to gross errors. Robust PCA (RPCA) has been proposed to remove the effect of sparse gross errors. For a given data matrix $M\in\mR^{m\times n}$, RPCA seeks to decompose it \sa{into} two parts $M:=\sa{L^o+S^o}$ \sa{where} $\sa{L^o}$ \sa{is} a low-rank matrix and $\sa{S^o}$ \sa{is} a sparse matrix. That is, RPCA assumes that $M$ is a superposition of $\sa{L^o}$ and $\sa{S^o}$. As a result, the gross errors will be captured by the sparse matrix $S^o$ so that the low-rank matrix $L^o$ can still approximate $M$ well. RPCA does not only provide a low-dimensional approximation which is robust to outliers, it also finds vast applications in a variety of real applications such as computer vision \cite{Candes-Li-Ma-Wright-RPCA}, \sa{image alignment} \cite{Peng-rasl-2012}, subspace recovery \cite{Liu-robust-subspace-recovery-2013}, clustering \cite{Shahid-iccv} and so on.

\sa{Mathematically, \cite{Candes-Li-Ma-Wright-RPCA,Chandrasekaran-RPCA-2011,NIPS2009_3704,NIPS2010_4005} investigated the conditions on the low-rank and sparse components $({L^o},{S^o})$ so that the inverse problem of recovering unknown $({L^o},{S^o})$ given $M$ is well defined. One particular formulation of RPCA can be stated as follows:}
\BE\label{rpca-rank-L0}
\min_{L,S\in\mR^{m\times n}} \rank(L) + \rho\|S\|_0 \ \st\ L+S = M,
\EE
where $\|S\|_0$ is called the $\ell_0$-norm\footnote{Technically, it is not a norm \sa{because it is not homogeneous;} but, we still call it a norm following the convention.} of $S$ and counts the number of nonzero entries of $S$, and $\rho>0$ is a
\sa{tradeoff} parameter. It is known that \eqref{rpca-rank-L0} is NP-hard and thus numerically intractable.
{Later,} it was shown in \cite{Candes-Li-Ma-Wright-RPCA,Chandrasekaran-RPCA-2011,NIPS2009_3704} that under certain conditions, \eqref{rpca-rank-L0} is equivalent to the following convex program with high probability:
\BE\label{rpca-nuclear-L1}
\min_{L,S\in\mR^{m\times n}} \|L\|_* + \rho\|S\|_1 \ \st\ L+S = M,
\EE
where $\|L\|_*$ is called the nuclear norm of $L$ and equal to the sum of the singular values of $L$, and $\|S\|_1:=\sum_{ij}|S_{ij}|$ is called the $\ell_1$ norm of $S$. \sa{The optimization problem in}~\eqref{rpca-nuclear-L1} is called robust principal component pursuit (RPCP), and it can be reformulated as a semidefinite program (SDP)~\sa{\cite{recht2010guaranteed}} and solved by \sa{an} interior point method \sa{for SDPs}. However, RPCA problems arising
{in} practice are usually of very large scale, and interior point methods do not scale well for these problems. More efficient algorithms that solve \eqref{rpca-nuclear-L1} and its variants {by exploiting the structure} in these problems were studied extensively in the literature. One variant of \eqref{rpca-nuclear-L1}
\sa{deals with an additional dense noise component. In particular, when}
\sa{$M$ contains also a} dense noise component {$N^o$ such that $\norm{N^o}_F\leq \sigma$ for some noise level $\sigma>0$}, i.e., $M=L^o+S^o+N^o$, \sa{instead of RPCP formulation} in \eqref{rpca-nuclear-L1},
the following so-called stable PCP (SPCP) problem \sa{is solved:}
\BE\label{rpca-nuclear-L1-noise}
\min_{L,S\in\mR^{m\times n}} \|L\|_* + \rho\|S\|_1 \ \st\ \norm{L+S-M}_F\leq\sigma.
\EE

{It is proved in \cite{Zhou-stable-pca-2010} that, {under certain conditions on $M$,} solving \eqref{rpca-nuclear-L1-noise} gives a stable estimate of $L^o$ and $S^o$ with high probability in the sense that $\|\hat{L}-L^o\|_F^2 + \|\hat{S}-S^o\|_F^2 \leq {\cO(mn\sigma^2)}$

where $(\hat{L},\hat{S})$ denotes the optimal solution to \eqref{rpca-nuclear-L1-noise}.}
\sa{Since \eqref{rpca-nuclear-L1-noise} satisfies the Slater's condition, it is equivalent to the following unconstrained problem for an appropriately chosen penalty parameter {$\mu>0$} depending on $\sigma$:
\BE\label{rpca-nuclear-L1-noise-reform-unconstrained}
\min_{L,S\in\mR^{m\times n}} \|L\|_* + \rho\|S\|_1 + \frac{{\mu}}{2} \|L+S-M\|_F^2.
\EE}

Note that if $M$ is only partially observed, that is, if we only have observations on $M_{ij}$ for some indices $(i,j)$ \sa{from} a subset $\Omega$, then \eqref{rpca-nuclear-L1}, \eqref{rpca-nuclear-L1-noise} and \eqref{rpca-nuclear-L1-noise-reform-unconstrained} can be respectively \sa{reformulated as}
\BE\label{rpca-nuclear-L1-partial}
\min_{L,S\in\mR^{m\times n}} \|L\|_* + \rho\|S\|_1 \ \st\ \sa{P_\Omega(L+S-M) = \mathbf{0},}
\EE
\BE\label{rpca-nuclear-L1-noise-partial}
\min_{L,S\in\mR^{m\times n}} \|L\|_* + \rho\|S\|_1 \ \st\ \sa{\norm{P_\Omega(L+S-M)}_F \leq \sigma,}
\EE
\BE\label{rpca-nuclear-L1-noise-reform-unconstrained-partial}
\min_{L,S\in\mR^{m\times n}} \|L\|_* + \rho\|S\|_1 + \frac{{\mu}}{2} \|P_\Omega(L+S-M)\|_F^2,
\EE
where the operator $P_\Omega: \mR^{m\times n}\rightarrow \mR^{m\times n}$ is defined as $[P_\Omega(M)]_{ij} = M_{ij}$, if $(i,j)\in\Omega$, and $[P_\Omega(M)]_{ij}=0$ otherwise. Most algorithms we discuss in this {paper} for solving \eqref{rpca-nuclear-L1}, \eqref{rpca-nuclear-L1-noise} and \eqref{rpca-nuclear-L1-noise-reform-unconstrained} can be used to solve \eqref{rpca-nuclear-L1-partial}, \eqref{rpca-nuclear-L1-noise-partial} and \eqref{rpca-nuclear-L1-noise-reform-unconstrained-partial} directly or with very little modification. For brevity, we will only describe algorithms for solving \eqref{rpca-nuclear-L1}, \eqref{rpca-nuclear-L1-noise} and \eqref{rpca-nuclear-L1-noise-reform-unconstrained} in this {paper}.

{For the sake of completeness,} \sa{we here briefly describe the results in \cite{Candes-Li-Ma-Wright-RPCA,Zhou-stable-pca-2010}. Assume that given data matrix $M\in\reals^{m\times n}$ is a superposition of the \emph{unknown} component matrices $L^o$, $S^o$, and $N^o$, i.e., $M=L^o+S^o+N^o$, such that $L^o$ is low-rank ($r:=\rank(L^0)\ll\min\{m,n\}$), $S^o$ is sparse ($s:=\norm{S^o}_0\ll mn$), and $\norm{N^o}\leq\sigma$ for some $\sigma>0$. Robust/stable PCA is an inverse problem with the objective of recovering $L^o$ and $S^o$ from the data matrix $M$. Clearly, if $S^o$ is low-rank and/or $L^o$ is sparse, the recovery is hopeless. To avoid these pathological instances, \cite{Candes-Li-Ma-Wright-RPCA,Zhou-stable-pca-2010} consider $(L^o,S^o)$ pairs coming from a particular class satisfying some incoherence and randomness conditions. Suppose the singular value decomposition (SVD) of $L^o$ is given by $L^o = \sum_{i=1}^r \sigma_iu_iv_i^\top$, where

$U=[u_1,\ldots,u_r]$ and $V = [v_1,\ldots,v_r]$ are formed by its left- and right-singular vectors. The incoherence conditions assume that there exists a parameter {$\delta>0$} such that
\begin{subequations}\label{incoherence-conditions}
\begin{align}
\max_i\|U^\top e_i\|^2\leq {\delta} r/m,\ & \max_i\|V^\top e_i\|^2 \leq {\delta} r/n, \label{incoherence-conditions-1} \\
\|UV^\top\|_\infty\leq\sqrt{\frac{{\delta} r}{mn}},\  & \label{incoherence-conditions-2}
\end{align}
\end{subequations}
where $\|Z\|_\infty := \max_{ij}|Z_{ij}|$ and $e_i$ denotes the $i$-th unit vector. For $S^o$, it is assumed that the set of indices for the nonzero entries is random and follows a uniform distribution among all the subsets of cardinality $s$. Roughly speaking, these conditions assure that the low-rank matrix $L^o$ is not sparse, and the sparse matrix $S^o$ is not low-rank.

Under these assumptions

it is shown in~\cite{Candes-Li-Ma-Wright-RPCA} that when $N^o=\mathbf{0}$, i.e., $\sigma=0$, solving the convex program \eqref{rpca-nuclear-L1} with $\rho=\frac{1}{\sqrt{\max\{m,n\}}}$ recovers the optimal solution of \eqref{rpca-rank-L0}, $(L^*,S^*)$, with high probability and $(L^*,S^*)=(L^o,S^o)$ provided that $L^o$ is sufficiently low-rank and $S^o$ is sufficiently sparse comparing to the matrix size -- see also~\cite{NIPS2009_3704}.} \sa{These results extend to the case where $M$ is partially observed; indeed, Cand\`{e}s et al.~\cite{Candes-Li-Ma-Wright-RPCA} show that solving~\eqref{rpca-nuclear-L1-partial} recovers $L^o$ under similar conditions.} Moreover, in \cite{Zhou-stable-pca-2010}, the authors showed that under certain conditions, solving the convex problem \eqref{rpca-nuclear-L1-noise} again with $\rho=\frac{1}{\sqrt{\max\{m,n\}}}$ \sa{generates a low-rank and sparse decomposition $(L^*,S^*)$ such that $\norm{L^*-L^o}_F^2+\norm{S^*-S^o}_F^2\leq Cmn{\sigma}^2$ for some {constant $C>0$ (independent of $m$, $n$ and $\sigma$)} with high probability -- note when $N=\mathbf{0}$, the recovery is exact with high probability.}

\red{Recently, there are works that further study statistical guarantees of different RPCA models. Zhang, Zhou and Liang \cite{Zhang-zhou-Liang} provide a refined analysis of RPCA which allows the support of the \sa{sparse} error matrix to be generated with non-uniform sampling, i.e., entries of the low-rank matrix are corrupted
with
different probabilities -- \sa{hence, one can model the scenario where some entries are more prone to corruption than the others}. A nonconvex model of RPCA is studied in \cite{Yi-Park-Chen-Caramanis-gradient} and a gradient descent method with proper initialization is shown to be able to reduce the computational complexity comparing with existing methods. Zhang and Yang \cite{Zhang-Yang-manifold-2017} consider a nonconvex optimization formulation with manifold constraint for RPCA. Two algorithms for manifold optimization are proposed in \cite{Zhang-Yang-manifold-2017} and it is shown that they can reduce the dependence on the \sa{condition} number of the underlying low-rank matrix theoretically. Netrapalli \etal \cite{Netrapalli-nonconvex-2014} consider another nonconvex formulation of RPCA and analyze the iteration complexity of the proposed alternating projection method.}

{There are also recent survey papers \cite{Bouwmans-survey-cviu-2014,Aybat-book-2016} that discuss algorithms for solving RPCA, but these papers mainly focus on its convex relaxations. In this {paper}, we aim to review existing algorithms for both convex and nonconvex relaxations/variants of RPCA models and point out a few possible future directions.}

{
\begin{table*}[t]
  {\centering
  \caption{{Convergence rates of different algorithms.}}\label{tab:convergence-rate}
   \begin{tabular}{llll}\thickhline
    \textbf{Algorithm}                                      & \textbf{Problem}                 & $\epsilon$-\textbf{optimality measure} & \textbf{Convergence rate}         \\\thickhline
    \multicolumn{4}{c}{\textbf{Algorithms for Convex Models}} \\\hline
    PGM \cite{Lin-WCAMSAP-2009}                    & \eqref{rpca-nuclear-L1-noise-reform-unconstrained} & objective value error  & $\cO(1/k)$            \\\hline
    APGM \cite{Lin-WCAMSAP-2009}                   & \eqref{rpca-nuclear-L1-noise-reform-unconstrained} & objective value error  & $\cO(1/k^2)$        \\\hline
    IALM \cite{Lin-Chen-Wu-Ma-iadm-2009}           & \eqref{rpca-nuclear-L1} & --- & convergence, no rate given \\\hline
    ADMM \cite{Yuan-Yang-pjo-2013}                 & \eqref{rpca-nuclear-L1} & --- & convergence, no rate given \\\hline
    ALM \cite{Goldfarb-Ma-Scheinberg-2013}         & \eqref{rpca-nuclear-L1-smooth-g} & objective value error & $\cO(1/k)$ \\\hline
    FALM \cite{Goldfarb-Ma-Scheinberg-2013}        & \eqref{rpca-nuclear-L1-smooth-g} & objective value error & $\cO(1/k^2)$ \\\hline
    ASALM \cite{Tao-Yuan-siopt-2011}               & \eqref{rpca-nuclear-L1-noise-N-rewrite} & --- & convergence unclear, no rate given \\\hline
    VASALM \cite{Tao-Yuan-siopt-2011}               & \eqref{rpca-nuclear-L1-noise-N-rewrite} & --- & convergence, no rate given \\\hline
    PSPG \cite{Aybat-Ma-Goldfarb-2014}             & \eqref{rpca-nuclear-L1-noise}    & objective value error & $\cO(1/k)$ \\\hline
    ADMIP \cite{Aybat-Iyengar-2015}                & \eqref{rpca-nuclear-L1-noise}    & {objective value error} & {$\cO(1/k)$} \\\hline
    Quasi-Newton method {(fastRPCA)}~\cite{P1C2-Aravkin14_1P}       & \eqref{eq:aravkin_model}         & --- & {convergence,} no rate given \\\hline
    3-block ADMM \cite{Lin-Ma-Zhang-2015-free-gamma} & \eqref{lin-ma-zhang-spcp}      & --- & convergence, no rate given \\\hline
    Frank-Wolfe \cite{Mu-FW-2016}                  & \eqref{FW-RPCP}                  & objective value error & $\cO(1/k)$ \\\hline
        \multicolumn{4}{c}{\textbf{Algorithms for Nonconvex Models}} \\\hline
    GoDec \cite{GoDec-icml-2011} & \eqref{zhou-tao} & ---  & local convergence, no rate given \\\hline
    GreBsmo \cite{pmlr-v31-zhou13b} & \eqref{eq:biliteral-model} & --- & convergence unclear, no rate given \\\hline
    Alternating Minimization (R2PCP) \cite{hintermuller2015robust} & \eqref{eq:ams-model} & --- & local convergence, no rate given \\\hline
    Gradient Descent (GD) \cite{Yi-Park-Chen-Caramanis-gradient} & $\approx$ \eqref{eq:biliteral-model} & --- & linear convergence with proper initialization \\ &&& and incoherence assumption \\\hline
    Alternating Minimization \cite{Gu-Wang-Liu-AISTATS-2016} & \eqref{gu-wang-liu} & --- & local convergence with proper initialization \\ &&& and incoherence and RIP assumptions \\\hline
    Stochastic alg. \cite{Xu-OR-PCA-2013} & \eqref{OR-PCA} & --- & convergence if the iterates are always full rank matrices, \\ &&& no rate given \\\hline
    LMafit \cite{Shen-Wen-Zhang-Lmafit} & \eqref{Lmafit-problem} & --- & convergence if difference between two consecutive iterates \\ &&& tends to zero, no rate given \\\hline
    Conditional Gradient \cite{Jiang-Lin-Ma-Zhang-nonconvex-2016} & \eqref{nonconvex-rpca-cg-constrained} & perturbed KKT & $\cO(1/\sqrt{k})$ \\\hline
    ADMM \cite{Jiang-Lin-Ma-Zhang-nonconvex-2016} & \eqref{nonconvex-rpca-admm} & perturbed KKT & $\cO(1/\sqrt{k})$ \\\hline
    Proximal BCD \cite{Jiang-Lin-Ma-Zhang-nonconvex-2016} & \eqref{nonconvex-rpca-bcd} & perturbed KKT & $\cO(1/\sqrt{k})$ \\\thickhline
    \end{tabular}
    }

    \rule{0in}{1.2em}$^\dag$\scriptsize {Note: Some of these algorithms solve different problems and the $\epsilon$-optimality measures are also different, so the convergence rates are not directly comparable with each other. \cite{Yi-Park-Chen-Caramanis-gradient} has no explicit optimization formulation but the objective is similar to \eqref{eq:biliteral-model}. Moreover, global convergence is usually guaranteed for convex solvers, but only local convergence is usually guaranteed for nonconvex solvers, unless certain very strong assumptions are made.}
\end{table*}
}

\section{Algorithms for Convex Relaxations/Variants of RPCA}\label{sec:convex}

The
{earliest first-order methods for} solving the convex RPCP problem
\sa{are given in}~\cite{Lin-WCAMSAP-2009,Lin-Chen-Wu-Ma-iadm-2009}. In \cite{Lin-WCAMSAP-2009}, the authors
proposed an accelerated proximal gradient method (APGM) \cite{Nesterov-07,Beck-Teboulle-2009,Tseng-2008} for solving \eqref{rpca-nuclear-L1-noise-reform-unconstrained} in which each iteration involves computing the proximal mappings of the nuclear norm $\|L\|_*$ and the $\ell_1$ norm $\|S\|_1$. In particular, the non-accelerated proximal gradient method (PGM) for solving \eqref{rpca-nuclear-L1-noise-reform-unconstrained} simply updates $L$ and $S$ as
\BE\label{pgm}
\BA{ll}
G^k     & := {\mu}(L^k+S^k-M) \\
L^{k+1} & := \argmin_L \ \|L\|_* + \frac{1}{2\tau}\|L-(L^k-\tau G^k)\|_F^2 \\
S^{k+1} & := \argmin_S \ \rho\|S\|_1 + \frac{1}{2\tau}\|S-(S^k-\tau G^k)\|_F^2.
\EA
\EE
Note that $G^k$ is the gradient of the quadratic penalty function in \eqref{rpca-nuclear-L1-noise-reform-unconstrained} and $\tau>0$ denotes a step size. The two subproblems in \eqref{pgm} both admit easy closed-form optimal solutions. Specifically, the solution of the $L$-subproblem corresponds to the proximal mapping of the nuclear norm, which is given by
\[L^{k+1} = \mathcal{S}_\tau(L^k-\tau G^k),\]
where the matrix shrinkage operation is defined as
\BE\label{mat-shrink}\mathcal{S}_\nu(Z) = U\diag((\sigma-\nu)_+)V^\top,\EE
where $Z = U\diag(\sigma)V^\top$ is the SVD of $Z$, and $z_+ := \max(0,z)$. The solution of the $S$-subproblem corresponds to the proximal mapping of the $\ell_1$ norm, which is given by
\[S^{k+1} = s_{\rho\tau}(S^k-\tau G^k),\]
where the vector shrinkage operation is defined as
\BE\label{vec-shrink}[s_\nu(Z)]_{ij} = \sign(Z_{ij})\circ\max\{0,|Z_{ij}|-\nu\},\EE
where $\sign(a)$ denotes the sign of $a$, and $\circ$ denotes the Hadamard product.

APGM incorporates Nesterov's acceleration technique and updates the variables as follows starting with $t_{-1}=t_0=1$.
\BE\label{apgm}
\BA{ll}
\bar{L}^k & := L^k + \frac{t_{k-1}-1}{t_k}(L^k-L^{k-1}) \\
\bar{S}^k & := S^k + \frac{t_{k-1}-1}{t_k}(S^k-S^{k-1}) \\[0.75mm]
\bar{G}^k & := {\mu}(\bar{L}^k+\bar{S}^k-M) \\
L^{k+1} & := \argmin_L \ \|L\|_* + \frac{1}{2\tau}\|L-(\bar{L}^k-\tau \bar{G}^k)\|_F^2 \\
S^{k+1} & := \argmin_S \ \rho\|S\|_1 + \frac{1}{2\tau}\|S-(\bar{S}^k-\tau \bar{G}^k)\|_F^2 \\
t_{k+1} & := (1+\sqrt{1+4t_k^2})/2.
\EA
\EE
Results in \cite{Nesterov-07,Beck-Teboulle-2009,Tseng-2008} show that the proximal gradient method \eqref{pgm} and the accelerated proximal gradient method \eqref{apgm} find an $\epsilon$-optimal solution to \eqref{rpca-nuclear-L1-noise-reform-unconstrained} in no more than $O(1/\epsilon)$ and $O(1/\sqrt{\epsilon})$ iterations, respectively.

\sa{When there is no noise, i.e., $\sigma=0$, the problem of interest is \eqref{rpca-nuclear-L1}. The drawback of the above approach for solving the unconstrained version in~\eqref{rpca-nuclear-L1-noise-reform-unconstrained} is that} \eqref{rpca-nuclear-L1-noise-reform-unconstrained} is equivalent to \eqref{rpca-nuclear-L1} only when ${\mu}\rightarrow +\infty$. Therefore, for any fixed ${\mu}>0$, there is always a residual term which does not go to zero. To remedy this, the same group of authors \cite{Lin-Chen-Wu-Ma-iadm-2009} considered the augmented Lagrangian method (ALM) for solving \eqref{rpca-nuclear-L1}. By associating a Lagrange multiplier $\Lambda$ to the linear equality constraint, the augmented Lagrangian function of \eqref{rpca-nuclear-L1} can be written as
\[\cL_\beta(L,S;\Lambda) := \|L\|_* + \rho\|S\|_1 - \langle \Lambda, L+S-M \rangle + \frac{\beta}{2}\|L+S - M\|_F^2,\]
where $\beta>0$ is a penalty parameter. A typical iteration of ALM {iterates the updates} as follows\sa{:}
\begin{subequations}
\label{rpca-nuclear-L1-ALM}
\begin{align}
(L^{k+1},S^{k+1}) & := \argmin_{L,S} \cL_\beta(L,S;\Lambda^k) \label{eq:ALM-subproblem}\\
\Lambda^{k+1} & := \Lambda^k - \beta(L^{k+1}+S^{k+1}-M).
\end{align}
\end{subequations}
Note that the first step in \eqref{rpca-nuclear-L1-ALM} requires to minimize the augmented Lagrangian function with respect to $L$ and $S$ simultaneously, which usually is \sa{computationally very expensive and almost as hard as}
solving the original problem \sa{in}~{\eqref{rpca-nuclear-L1}}. In \cite{Lin-Chen-Wu-Ma-iadm-2009}, the authors proposed both exact and inexact versions of ALM, where the former one solves the subproblems (almost) exactly and the latter one solves the subproblems inexactly
\sa{according to a particular subproblem termination criterion.}
Both the exact ALM and inexact ALM (IALM) employ some iterative algorithm for minimizing the augmented Lagrangian function until certain \sa{overall} stopping criterion is met, which may require many iterations and thus time consuming.
Around the same time when \cite{Lin-WCAMSAP-2009,Lin-Chen-Wu-Ma-iadm-2009} appeared, the alternating direction method of multipliers (ADMM) was revisited and found very successful in solving signal processing and image processing problems \cite{Combettes-2007,Goldstein-Osher-split-Bregman-2009,Yang-Zhang-Yin-ADMM-2010,Yang-Zhang-admm-2011}. It was then found that RPCP \sa{in} \eqref{rpca-nuclear-L1} can be nicely solved by ADMM due to its special separable structure \cite{Yuan-Yang-pjo-2013,Goldfarb-Ma-Scheinberg-2013}.

\sa{The ADMM iterations for solving \eqref{rpca-nuclear-L1} take the following form:}
\BEA\label{rpca-nuclear-L1-ADMM}
\BA{ll}
L^{k+1} & := \argmin_{L} \cL_\beta(L,S^k;\Lambda^k) \\
S^{k+1} & := \argmin_{S} \cL_\beta(L^{k+1},S;\Lambda^k) \\
\Lambda^{k+1} & := \Lambda^k - \beta(L^{k+1}+S^{k+1}-M).
\EA
\EEA
Comparing to ALM \sa{in} \eqref{rpca-nuclear-L1-ALM}, it is noted that ADMM splits the subproblem in \sa{\eqref{eq:ALM-subproblem} into} two \sa{smaller} subproblems that correspond to computing proximal mappings of $\|L\|_*$ and $\|S\|_1$, respectively. The ADMM \eqref{rpca-nuclear-L1-ADMM} is known as two-block ADMM as there are two block variables $L$ and $S$ and hence two subproblems \sa{are solved} in each iteration of the algorithm.
It is now widely known that the two-block ADMM is a special case of the so-called Douglas-Rachford operator splitting method \cite{Gabay-Mercier-1976,Gabay-83,Fortin-Glowinski-1983,Lions-Mercier-79} applied to the dual problem, and

the two-block ADMM for solving convex problems globally converges for any penalty parameter $\beta>0$ \cite{Eckstein-Bertsekas-1992} and converges with a sublinear rate $O(1/k)$ (see, e.g., \cite{He-Yuan-rate-ADM-2012,Monteiro-Svaiter-2010a,Lin-Ma-Zhang-convergence-2014}).

The alternating linearization method (ALM) proposed by Goldfarb, Ma and Scheinberg \cite{Goldfarb-Ma-Scheinberg-2013} is shown to be equivalent to a symmetric version of ADMM \eqref{rpca-nuclear-L1-ADMM} with either $\|L\|_*$ or $\|S\|_1$ \sa{replaced with some suitable smooth approximation}. \sa{For instance, given $\nu>0$, define $g_\nu:\reals^{m\times n}\rightarrow\reals$ such that
\begin{align}
\label{eq:smooth-g}
g_{\nu}(S)=\max_{Z\in\reals^{m\times n}}\{\fprod{S,Z}-\frac{\nu}{2}\norm{Z}_F^2:\ \norm{Z}_{\infty}\leq \rho\},
\end{align}
and let $g(S)=\rho \norm{S}_1$. Clearly, $g_\nu\rightarrow g$ uniformly as $\nu\searrow 0$. Moreover, $g_\nu$ is a differentiable convex function such that $\grad g_\nu$ is Lipschitz continuous. Indeed, given $S\in\reals^{m\times n}$, let $Z_\nu(S)$ be the maximizer for \eqref{eq:smooth-g}, which in closed form can be written as $Z_\nu(S)={\sign(S)\circ}\max\{\frac{1}{\nu}|S|,~\rho \mathbf{1}_{m\times n}\}$, and $\grad g_\nu(S)=Z_\nu(S)$ is Lipschitz {continuous} with constant $C_{g_{\nu}}=\frac{1}{\nu}$. Similarly, given $\mu>0$, define $f_\mu:\reals^{m\times n}\rightarrow\reals$ such that
\begin{align}
\label{eq:smooth-f}
f_{{\mu}}(L)=\max_{W\in\reals^{m\times n}}\{\fprod{L,W}-\frac{\mu}{2}\norm{W}_F^2:\ \norm{W} \leq 1\},
\end{align}
where $\norm{\cdot}$ denotes the spectral norm, and let $f(L)=\lambda \norm{L}_*$. Clearly, $f_\mu\rightarrow f$ uniformly as $\mu\searrow 0$. Moreover, $f_\mu$ is a differentiable convex function such that $\grad f_\mu$ is Lipschitz continuous. Indeed, given $L\in\reals^{m\times n}$, let $W_\mu(L)$ be {the maximizer for \eqref{eq:smooth-f},} which in closed form can be written as \magenta{$W_\mu(L)=U\diag\big(\max\big\{\frac{\sigma}{\mu}-\mathbf{1},~\mathbf{0}\big\}\big)V^\top$,} where $L=U \diag(\sigma)V^\top$ is the singular value decomposition of $L$ with $\sigma\in\reals^r_{++}$ denoting the vector of singular values; moreover, $\grad f_\mu(S)=W_\mu(L)$ is Lipschitz {continuous} with constant $C_{f_{\mu}}=\frac{1}{\mu}$.}

The alternating linearization method in \cite{Goldfarb-Ma-Scheinberg-2013} can be applied to solve the following problem, which is a smoothed version of \eqref{rpca-nuclear-L1}:
\BE\label{rpca-nuclear-L1-smooth-g} \min \ \|L\|_* + g_{\nu}(S), \ \st, \ \ L+S = M.\EE
Denote the augmented Lagrangian function of \eqref{rpca-nuclear-L1-smooth-g} as $\tilde{\cL}_\beta(L,S;\Lambda)$, the alternating linearization method in  \cite{Goldfarb-Ma-Scheinberg-2013} {iterates the updates} as follows:
\BE\label{alternating-linearization-method}
\BA{ll}
L^{k+1} & := \argmin_{L} \tilde{\cL}_\beta(L,S^k;\Lambda^k), \\
\Lambda^{k+\half} & := \Lambda^k - \beta(L^{k+1}+S^{k}-M), \\
S^{k+1} & := \argmin_{S} \tilde{\cL}_\beta(L^{k+1},S;\Lambda^{k+\half}), \\
\Lambda^{k+1} & := \Lambda^{k+\half} - \beta(L^{k+1}+S^{k+1}-M).
\EA
\EE
The authors in \cite{Goldfarb-Ma-Scheinberg-2013} proved that this method has a sublinear convergence rate $O(1/k)$. They also proposed an accelerated version of \eqref{alternating-linearization-method} (FALM) in \cite{Goldfarb-Ma-Scheinberg-2013} by adopting Nesterov's acceleration technique, and proved that the accelerated alternating linearization method has a better sublinear convergence rate $O(1/k^2)$.

Based on the success of two-block ADMM for solving \eqref{rpca-nuclear-L1}, it is then very natural to apply ADMM to solve SPCP in \eqref{rpca-nuclear-L1-noise}. {To do so, one has to introduce a new variable $N$, and rewrite \eqref{rpca-nuclear-L1-noise} equivalently as
\BE\label{rpca-nuclear-L1-noise-N-rewrite}
\BA{ll}
\min_{L,S,N\in\mR^{m\times n}} & \|L\|_* + \rho\|S\|_1 \\ \st, & L+S+N=M, \norm{N}_F\leq\sigma.
\EA
\EE
The ADMM for solving \eqref{rpca-nuclear-L1-noise-N-rewrite} {iterates} as follows with three block variables:
\BEA\label{spcp-ADMM}
\BA{ll}
L^{k+1} & := \argmin_{L} \cL_\beta(L,S^k,N^k;\Lambda^k) \\
S^{k+1} & := \argmin_{S} \cL_\beta(L^{k+1},S,N^k;\Lambda^k) \\
N^{k+1} & := \argmin_{N} \cL_\beta(L^{k+1},S^{k+1},N;\Lambda^k) \\
\Lambda^{k+1} & := \Lambda^k - \beta(L^{k+1}+S^{k+1}+N^{k+1}-M),
\EA
\EEA
where the augmented Lagrangian function for \eqref{rpca-nuclear-L1-noise-N-rewrite} is defined as}
\sa{
\[\BA{l}\cL_\beta(L,S,N;\Lambda) := \|L\|_* + \rho\|S\|_1 + \Ind(N\mid\|N\|_F\leq\sigma) \\ - \langle \Lambda, L+S+N-M \rangle + \frac{\beta}{2}\|L+S+N-M\|_F^2,\EA\]}
where $\Ind(N\mid\cN)$ denotes the indicator function of the set $\{N\in\cN\}$, i.e., $\Ind(N\mid\cN)=0$ if $N\in\cN$ and $\Ind(N\mid\cN)=\infty$ otherwise. Note that the three subproblems in \eqref{spcp-ADMM} all have closed-form solutions. In particular, the $L$-subproblem corresponds to the proximal mapping of $\|L\|_*$, the $S$-subproblem corresponds to proximal mapping of $\|S\|_1$, and the $N$-subproblem corresponds to projection onto the set $\{N\mid \|N\|_F\leq\sigma\}$. Similar idea was used in \cite{Peng-rasl-2012} for robust image alignments. In practice, this three-block ADMM usually works very well. However, it was later discovered that the ADMM with more than two block variables is not necessarily convergent in general \cite{Chen-admm-failure-2013}. \sa{Note that although {\eqref{rpca-nuclear-L1-noise-N-rewrite}} contains three block variables, it can be viewed as a two-block problem, if we group $S$ and $N$ as one (larger) block variable. One of the earliest methods for solving SPCP in \eqref{rpca-nuclear-L1-noise} and \eqref{rpca-nuclear-L1-noise-N-rewrite} is a three-block ADMM algorithm, ASALM, proposed by Tao and Yuan~\cite{Tao-Yuan-siopt-2011}, and although it does not have any convergence guarantees, it works well in practice; and slightly changing the update rule in ASALM leads to VASALM, of which iterate sequence converges to an optimal solution; but this comes at the cost of degradation in practical convergence speed when compared to ASALM -- {indeed, VASALM \cite{Tao-Yuan-siopt-2011} can be seen as a linearized version of two-block ADMM to solve \eqref{rpca-nuclear-L1-noise} with a convergence guarantee without any convergence rate result.} To remedy the shortcoming associated with the theoretical convergence of three-block ADMM, several other alternatives based on two-block ADMM were proposed~\cite{aybat2011fast,Aybat-Iyengar-2015,Aybat-Ma-Goldfarb-2014}.}

\sa{Aybat, Goldfarb and Ma \cite{Aybat-Ma-Goldfarb-2014} proposed an accelerated proximal gradient method, \sa{PSPG,} for solving SPCP in~\eqref{rpca-nuclear-L1-noise}.
First, \eqref{rpca-nuclear-L1-noise} is reformulated with a partially smooth objective. In particular, the nuclear norm is smoothed according to \eqref{eq:smooth-f}:

\begin{align}
\label{rpca-nuclear-L1-noise-rewrite-smooth}
\min_{L,S\in\mR^{m\times n}} f_\mu(L) + \rho\|S\|_{1}\ \st\ (L,S)\in\chi,\\
\chi:=\{(L,S)\mid\ \|L+S - M\|_F \leq \sigma\},
\end{align}
where $\mu>0$ is a given smoothing parameter. An accelerated proximal gradient method such as~\cite{Beck-Teboulle-2009,Nesterov-07}} can be applied to solve \eqref{rpca-nuclear-L1-noise-rewrite-smooth}, because it was shown in \cite{Aybat-Ma-Goldfarb-2014} that the following subproblem is easy to solve:
\begin{align}
\label{eq:complicated-prox}
\min_{L,S} \|S\|_1+\frac{1}{2\xi}\|L-\tilde{L}\|_F^2\ \st\ (L,S)\in\chi,
\end{align}
where $\xi>0$ denotes a step size of the proximal gradient step and $\tilde{L}$ denotes some known matrix. \sa{This operation requires one sorting which has $\cO(mn\log(mn))$ complexity.}

\sa{For any $\epsilon>0$, setting $\mu=\Omega(\epsilon)$, PSPG proposed in~\cite{Aybat-Ma-Goldfarb-2014} can compute an $\epsilon$-optimal solution to \eqref{rpca-nuclear-L1-noise} within $\cO(1/\epsilon)$ iterations, and its computational complexity per iteration is comparable to the work per iteration required by ASALM and VASALM, which is mainly determined by an SVD computation. On the other hand, it is also important to emphasize that PSPG iterate sequences do not converge to an optimal solution to the SPCP problem in \eqref{rpca-nuclear-L1-noise}. In particular, since within PSPG the smoothing parameter $\mu$ is fixed, depending on the approximation parameter $\epsilon$ for solving \eqref{rpca-nuclear-L1-noise-rewrite-smooth}, further iterations after reaching an $\epsilon$-optimal solution in $\cO(1/\epsilon)$ iterations do not necessarily improve the solution quality.}

In~\cite{aybat2011fast,Aybat-Iyengar-2015}, the \emph{variable penalty} ADMM algorithm ADMIP (Alternating Direction Method with Increasing Penalty) is proposed to solve the following equivalent formulation for \eqref{rpca-nuclear-L1-noise} using the variable splitting trick:
\BE\label{rpca-nuclear-L1-noise-equiv}
\min_{\hat{L},L,S\in\mR^{m\times n}} \|\hat{L}\|_* + \rho\|S\|_1 \ \st\ (L,S)\in\chi,\quad L=\hat{L}.
\EE
The augmented Lagrangian function of \eqref{rpca-nuclear-L1-noise-equiv} can be written as
\[\cL_\beta(\hat{L},L,S;\Lambda) := \|\hat{L}\|_* + \rho\|S\|_1 - \langle \Lambda, \hat{L}-L \rangle + \frac{\beta}{2}\|\hat{L}-L\|_F^2.\]
Given a nondecreasing penalty parameter sequence $\{\beta^k\}_{k\in\integers_+}$, ADMIP {updates the variables} as follows\sa{:}
\begin{subequations}
\label{rpca-nuclear-L1-ADMIP}
\begin{align}
\hat{L}^{k+1} & := \argmin_{\hat{L}} \cL_{\beta^k}(\hat{L},L^k,S^k;\Lambda^k) \label{eq:ADMIP-subproblem-1}\\
(L^{k+1},S^{k+1}) & := \argmin_{(L,S)\in\chi} \cL_{\beta^k}(\hat{L}^{k+1},L,S;\Lambda^k) \label{eq:ADMIP-subproblem-2}\\
\Lambda^{k+1} & := \Lambda^k - \beta^k(\hat{L}^{k+1}-L^{k+1}).
\end{align}
\end{subequations}
The step in~\eqref{eq:ADMIP-subproblem-1} requires computing a soft thresholding on the singular values of an $m\times n$ matrix and the step in~\eqref{eq:ADMIP-subproblem-2} requires an operation given in \eqref{eq:complicated-prox}.

Under mild conditions on the penalty parameter sequence, Aybat and Iyengar show that the primal-dual ADMIP iterate sequence converges to an optimal primal-dual solution to the SPCP problem in \eqref{rpca-nuclear-L1-noise-equiv} -- hence, $\{(L^k,S^k)\}_{k\in\integers_+}$ converges to an optimal solution to \eqref{rpca-nuclear-L1-noise}, and when constant penalty parameter is used as a special case, it can compute an $\epsilon$-optimal solution within $\cO(1/\epsilon)$ iterations, of which complexity is determined by an SVD. In particular, one needs the penalty parameter sequence $\{\beta^k\}_{k \in\integers_+}$ to be non-decreasing and to satisfy $\sum_{k}(\beta^k)^{-1}=+\infty$. The main advantages of adopting an increasing sequence of penalties are as follows:
  \begin{enumerate}
  \item The algorithm is robust in the sense that there is no need to search for problem data dependent $\beta^*$ that works well in practice.
  \item The algorithm is likely to achieve primal feasibility faster.
  \item The complexity of initial (transient) iterations can be controlled through controlling $\{\beta^k\}$.

    The main computational bottleneck in ADMIP is the SVD computation in \eqref{eq:ADMIP-subproblem-1}. Since the optimal $L^*$ is of low-rank, and $L_k\rightarrow L^*$, eventually the SVD computations are likely to be very efficient. However, since the initial iterates in the transient phase of the algorithm may have large rank, the complexity of the SVD in the initial iterations
    can be quite large. To compute the solution to the subproblem in~\eqref{eq:ADMIP-subproblem-1}, one does not need to compute singular values smaller than $1/\beta^k$; hence, initializing ADMIP with a \emph{small} $\beta^0>0$ will significantly decrease the complexity of initial iterations through employing \emph{partial} SVD computations, e.g., Lanczos-based methods such as PROPACK~\cite{larsen1998lanczos}.
  \end{enumerate}
In~\cite{Aybat-Iyengar-2015}, Aybat and Iyengar compared ADMIP against ASALM on both randomly generated synthetic problems and surveillance video foreground extraction problems. According to numerical results reported in~\cite{Aybat-Iyengar-2015}, on the synthetic problems ASALM requires about \emph{twice} as many iterations for convergence, while the total runtime for ASALM is considerably larger.

\sa{Aravkin et al.~\cite{P1C2-Aravkin14_1P} proposed solving
\begin{equation}
\label{eq:aravkin_model}
\min\psi(L,S)\quad \hbox{ s.t. }\quad \phi(\cA(L,S)-M) \leq \sigma,
\end{equation}
where $\cA:\reals^{m\times n}\times\reals^{m\times n}\rightarrow\reals^{m\times n}$ is a linear operator, $\phi:\reals^{m\times n}\rightarrow\reals$ is a smooth convex loss-function, and $\psi$ can be set to either one of the following functions:
\begin{align*}
\psi_{\rm sum}(L,S)&:=\norm{L}_* + \rho\norm{S}_1,\\
\psi_{\rm max}(L,S)&:=\max\{\norm{L}_*, \rho_{\max}\norm{S}_1\},
\end{align*}
$\rho,\rho_{\max}>0$ are some given function parameters. Note that setting $\psi=\psi_{\rm sum}$, $\rho(.)=\norm{.}^2_F$, and $\cA(L,S)=\pi_{\Omega}{L+S}$ in \eqref{eq:aravkin_model}, one obtains the SPCP problem in \eqref{rpca-nuclear-L1-noise-partial}. This approach offers advantages over the original SPCP formulation in terms of practical parameter selection. The authors make a case that although setting $\rho=\frac{1}{\sqrt{\max\{m,n\}}}$ in \eqref{rpca-nuclear-L1-noise} has theoretical justification as briefly discussed in the introduction section, many practical problems may violate the underlying assumptions in \eqref{incoherence-conditions}; in those cases one needs to tune $\rho$ via cross validation, and selecting $\rho_{\max}$ in $\psi_{\max}$ might be easier than selecting $\rho$ in $\psi_{\rm sum}$. Instead of solving \eqref{eq:aravkin_model} directly, a convex variational framework, accelerated with a ``quasi-Newton" method, is proposed. In particular, Newton's method is used to find a root of the value function:
\BE\label{eq:aravkin_value_function}
\BA{ll}
\upsilon(\tau):= & \min_{L,S\in\reals^{m\times n}}\phi(\cA(L,S)-A)-\sigma \\
                 & \st, \psi(L,S)\leq\tau,
\EA
\EE
i.e., given $\sigma>0$ compute $\tau^*=\tau(\sigma)$ such that $\upsilon(\tau^*)=0$. According to results in~\cite{P1C2-Aravkin13_1J}, if the constraint in \eqref{eq:aravkin_model} is tight at an optimal solution, then there exists $\tau^*=\tau(\sigma)$ such that $\upsilon(\tau^*)=0$ and the corresponding optimal solution to \eqref{eq:aravkin_value_function} is also optimal to \eqref{eq:aravkin_model}. Within Newton's method for root finding, to compute the next iterate $\tau^{k+1}$, one can compute the derivative of the value function at the current iterate $\tau^k$ as follows $\upsilon'(\tau_k)=-\psi^\circ(\cA^\top\grad\phi(\cA(L_k,S_k)-A))$, where $\psi^\circ$ denotes the polar gauge to $\psi$ and $(L_k,S_k)$ denotes the optimal solution to \eqref{eq:aravkin_value_function} at $\tau=\tau_k$ -- Aravkin et al. proposed a projected ``Quasi-Newton" method to solve \eqref{eq:aravkin_value_function}. According to numerical tests reported in \cite{P1C2-Aravkin13_1J}, QN-max, the quasi-Newton method running on \eqref{eq:aravkin_model} with $\psi=\psi_{\rm max}$ and $\phi(.)=\norm{.}_F^2$, is competitive with the state-of-the-art codes, ASALM~\cite{Tao-Yuan-siopt-2011}, PSPG~\cite{Aybat-Ma-Goldfarb-2014}, and ADMIP~\cite{Aybat-Iyengar-2015}.}

{
In a recent work \cite{Lin-Ma-Zhang-2015-free-gamma}, Lin, Ma and Zhang considered \sa{the penalty formulation of the SPCP problem, which is equivalent to solving \eqref{rpca-nuclear-L1-noise} for certain noise level $\sigma>0$:}
\BE\label{lin-ma-zhang-spcp}
\BA{ll}
\min & \|L\|_* + \rho\|S\|_1 + {\mu}\|\sa{N}\|_F^2 \\
\st  & L + S + \sa{N} = M,
\EA
\EE
where $\rho>0$ is the sparsity tradeoff parameter and $\mu>0$ is a suitable penalty parameter depending on the noise level $\sigma>0$. The authors showed that the following 3-block ADMM for solving \eqref{lin-ma-zhang-spcp} globally converges for any penalty parameter ${\beta}>0$.
{
\BE\label{lin-ma-zhang-spcp-admm}
\BA{ll}
L^{k+1} & := \argmin_L \ \cL_\beta(L,S^k,N^k;\Lambda^k) \\
S^{k+1} & := \argmin_S \ \cL_\beta(L^{k+1},S,N^k;\Lambda^k) \\
N^{k+1} & := \argmin_N \ \cL_\beta(L^{k+1},S^{k+1},N;\Lambda^k) \\
\Lambda^{k+1} & := \Lambda^k - \beta(L^{k+1} + S^{k+1} + N^{k+1} - M),
\EA
\EE}
where the augmented Lagrangian function is
\sa{
\[\BA{l}\cL_\beta(L,S,N;\Lambda) := \|L\|_* + \rho\|S\|_1 + \mu\|N\|_F^2 \\ - \langle \Lambda, L + S + N - M\rangle + \frac{{\beta}}{2}\|L + S + N - M\|_F^2.\EA\]}
Note that the three subproblems in \eqref{lin-ma-zhang-spcp-admm} are all easy to solve. Specifically, the $L$-subproblem corresponds to the proximal mapping of nuclear norm $\|L\|_*$, the $S$-subproblem corresponds to the proximal mapping of $\ell_1$ norm, and the $N$-subproblem admits a very easy analytical solution.
}

The Frank-Wolfe method (aka conditional gradient method) \cite{FW-1956} was revisited recently for solving large-scale machine learning problems \cite{Jaggi-icml-2013,thesisMJ}. RPCA is a representative example that is suitable for Frank-Wolfe method. Note that algorithms discussed above usually involve computing the proximal mapping of the nuclear norm, which is given by an SVD in \eqref{mat-shrink}. Computing full SVD for a large matrix in every iteration can be very time consuming. In contrast, the Frank-Wolfe method deals with nuclear norm in a much simpler manner, which only computes the largest singular value of a matrix in each iteration. The Frank-Wolfe method for solving RPCA was proposed by Mu \etal in \cite{Mu-FW-2016}. The authors in \cite{Mu-FW-2016} considered the penalized variant of RPCA \eqref{rpca-nuclear-L1-noise-reform-unconstrained}. However, this problem cannot be directly solved by the Frank-Wolfe method, because the Frank-Wolfe method requires a bounded constraint set. Therefore, the authors further reformulated \eqref{rpca-nuclear-L1-noise-reform-unconstrained} to the following problem for properly chosen constants $\lambda_L$, $\lambda_S$, $U_L$ and $U_S$:
\BE\label{FW-RPCP}
\BA{ll}
\min  & \half\|L+S-M\|_F^2 + \lambda_L t_L + \lambda_S t_S, \\
\st,                         & \|L\|_*\leq t_L \leq U_L,\ \|S\|_1\leq t_S \leq U_S.
\EA
\EE
The Frank-Wolfe method {iterates the updates} as follows:
\BE\label{FW-alg}
\BA{ll}
G^k & := L^k + S^k -M \\
\sa{(d_L^k,d_{t_L}^k)} & := \argmin_{\|L\|_*\leq{t_L}\leq U_L} \ \langle G^k,~L\rangle + \lambda_L{t_L}\\
\sa{(d_S^k,d_{t_S}^k)} & := \argmin_{\|S\|_1\leq{t_S}\leq U_S} \ \langle G^k,~S\rangle + \lambda_S{t_S} \\
\gamma^k &:= 2/(k+2)\\
L^{k+1} & := (1-\gamma^k)L^k + \gamma^k d_L^k  \\
t_L^{k+1} & := (1-\gamma^k)t_L^k + \gamma^k d_{t_L}^k \\
S^{k+1} & := (1-\gamma^k)S^k + \gamma^k d_S^k  \\
t_S^{k+1} & = (1-\gamma^k)t_S^k + \gamma^k d_{t_S}^k.
\EA
\EE
It was shown in \cite{Mu-FW-2016} that the two minimization subproblems in \eqref{FW-alg} are easy to solve. In particular, solving the subproblem for $(d_L,d_{t_L})$ requires only to compute the largest singular value and its corresponding singular vector of an $m\times n$ matrix. This is a big saving compared with computing the full SVD as required for computing the proximal mapping of the nuclear norm. As a result, Frank-Wolfe method has better per-iteration complexity than the proximal gradient method and ADMM algorithms discussed above, and thus may have better scalability for very large-scale problems.
\sa{On the other hand, as pointed out in~\cite{Mu-FW-2016}, one clear disadvantage of Frank-Wolfe method on \eqref{FW-RPCP} is that at every iteration only one entry of the sparse component is updated. This leads to very slow convergence in practice. Hence, Mu et al.~\cite{Mu-FW-2016} proposed combining Frank-Wolfe iterations with an additional proximal gradient step in $S$-block. In particular, they proposed after Frank-Wolfe iterate $(L_{k+1},\tilde{S}_{k+1})$ is computed, an extra proximal gradient step is computed and $S$-block is updated again. Moreover, the authors also showed that this hybrid method obtained by combining Frank-Wolfe and proximal gradient steps enjoys a sublinear convergence rate $O(1/k)$ similar to Frank-Wolfe method given in~\eqref{FW-alg}.}

As a special case of RPCA, one can consider that all columns of the low-rank matrix $L$ are identical. That is, the given matrix $M$ is a superposition of a special rank-one matrix $L$ and a sparse matrix $S$. This special RPCA finds many interesting applications in practice such as video processing \cite{Li-Ng-Yuan-NLAA-2015,Li-Wang-Hu-Cai-2014} and bioinformatics \cite{Ma-Sparcoc-2015}. For instance, in the background extraction of surveillance video, if the background is static, then the low-rank matrix $L$ that corresponds to the background should have identical columns. As a result, the background and foreground can be separated by solving the following convex program:
\BE\label{nuclear-norm-free}
\BA{ll}
\min_{x,S} & \|S\|_1 \\
\st & [x,x,\cdots,x] \circ \mathbf{E} + S = M,
\EA
\EE
where $[x,x,\cdots,x]$ denotes the $m\times n$ matrix with all columns being $x$, $\mathbf{E}$ denotes the $m\times n$ {matrix} with all ones. Note that the optimal $x$ of \eqref{nuclear-norm-free} corresponds to the static background for all frames and $S$ corresponds to the moving foreground. The advantage of \eqref{nuclear-norm-free} is that it does not involve nuclear norm. As a result, SVD can be avoided when designing algorithms for solving it which makes the resulting algorithms very efficient. Yang, Pong and Chen \cite{Yang-Pong-Chen-2017} adopted the similar idea and designed variants of ADMM algorithm for solving a more general model where the sparsity function of $S$ is allowed to be a nonconvex function. Convergence of the proposed ADMM was proved under the assumption of KL property \cite{Kurdyka-1998,Lojasiewicz-1963} being satisfied. {We will discuss these topics in more details in the next section.}


TITLE: Inequality for the $L^2$ norm of a derivative in a Dirichlet problem
QUESTION [2 upvotes]: I found some trouble in solving this problem:
given $f$ a continuous function in $[0,1]$ and $u\in C^2([0,1])$, such that $u''(x)=f(x)$ on $(0,1)$ and $u(0)=u(1)=0$, prove, $\forall\epsilon > 0$, the inequality
$||u'||_2^2 \leq \epsilon||u||_2^2 + \frac{1}{4\epsilon}||f||_2^2$
where $||.||_2$ is the usual norm on $L^2 ([0,1])$.
First, I rewrote the equation in this way:
$-u(x)u''(x)=-u(x)f(x)$
in order to have
$||u'||_2^2 = -\int_{0}^{1}u(x)f(x)dx$.
Then I tried to use a Green function to represent $u$, i.e.
$F(x,y)$=\begin{cases} -x(1-y) & \text{if $0\leq x\leq y\leq 1$} \\ -y(1-x) & \text{if $0\leq y\leq x\leq 1$} \end{cases}
in order to have $u(x)=\int_{0}^{1}F(x,y)f(y)dy$ but I really don't know where that $\epsilon$ should come from. I also thought I could use some identity approximation or mollifiers, but I got stuck.

REPLY [1 votes]: Hint:
$$ \epsilon  u^2 + \frac{1}{4\epsilon} f^2 + uf = \left( \sqrt{\epsilon} u + \frac{1}{2\sqrt{\epsilon}}f \right)^2 $$


TITLE: In sage: from $\mathbb{Z}/N\mathbb{Z}$ to $(\mathbb{Z}/N\mathbb{Z})^*$
QUESTION [3 upvotes]: In Sage G=IntegerModRing(11^5-1) gives me the cyclic group $\mathbb{Z}/161050\mathbb{Z}$. From here U=G.unit_group() gives me the group $(\mathbb{Z}/161050\mathbb{Z})^*\simeq\mathbb{Z}/20\mathbb{Z}\times\mathbb{Z}/3220\mathbb{Z}$ .
Question: how can I get in Sage the element in $\mathbb{Z}/20\mathbb{Z}\times\mathbb{Z}/3220\mathbb{Z}$ corresponding to $11\in\mathbb{Z}/161050\mathbb{Z}$?
There is a map from $(\mathbb{Z}/161050\mathbb{Z})^*$ to $\mathbb{Z}/161050\mathbb{Z}$: we can do U.inject_variables() to get variables f0 and f1 and then make G(f0) and G(f1) to get the corresponding elements of $\mathbb{Z}/161050\mathbb{Z}$. But U(G(11)) responds an error.
Of course the number here a only example and I search a general method.

REPLY [4 votes]: U(G(11).generalised_log()) does what you want. See the documentation for this method here.


TITLE: What is the mechanism behind bathtub whirl?
QUESTION [1 upvotes]: I understand that at the middle of the whirl can't be any water because that would make the velocity of water infinite.
But what is the mechanism, the force, the structure that removes water from the center?

REPLY [1 votes]: For low draining rates you can have water at the centre of the whirl. The free vortex model in which tangential velocity of the whirl goes as $1/r$ (where $r$ is the radial distance from the whirl's centre) is an idealization. In that model you get infinite velocity at the centre ($r=0$) because it is assumes that infinite amount of vorticity is concentrated at the centre of the whirl. But real fluids have viscosity which will not allow such infinite concentration and instead smear the vorticity over a finite radius. A better model for such whirls is a Rankine vortex which contains a core of radius $R$ which is in solid body rotation and therefore in this region the velocity goes as $r$, while outside this core the velocity goes as $1/r$. Thus the singularity is avoided by nature.


TITLE: Error estimation in Taylor Series
QUESTION [2 upvotes]: $$f(x) = \ln (1-x) $$
I should find third power series solution in point $0$. First, I find the derivatives:
$$ f'(x) = -\frac1{1-x} $$
$$ f''(x) = -\frac1{(1-x)^2} $$
$$ f'''(x) = -\frac2{(1-x)^3}$$
Now put them in the formula:
$$ f(x) \approx \ln 1 - \frac11(x-0) - \frac{1}{2!}(x-0)^2 - \frac{2}{3!}(x-0)^3 = \ln 1 - x - \frac{x^2}{2} - \frac{2x^3}{6}$$
I should now estimate it's error if $|x| \leq 0.5$
What exactly should I do after the formula part.

REPLY [2 votes]: Some observations:
Calculate the next few terms, and reduce the fractions to lowest terms. Do you see a pattern? Can you see why it's there?
If $x$ is negative, how do the signs behave? How can you use that to get an error bound really easily (when $x$ is negative)?
What do you know about the value of the Taylor remainder?


TITLE: What does it mean when two Groups are isomorphic?
QUESTION [38 upvotes]: I'm not asking for the formal definition I know it. An isomorphism is a bijective homomorphism. In my book it's indicated many times when two groups are isomorphic, and I don't understand what's the reason for that. What can we "do" when we know that 2 groups are isomorphic? What does it really mean when 2 groups are isomorphic?

REPLY [0 votes]: As I have commented below the answers Bamboo and Ethan Bolker, ismomorphism between two groups  does not mean "they are one and the same".  They are same only theoretically, i.e., algebraically, but not practically. Actually they could be differnet practically is the reason isomorphisms are useful!
Isomorphism only means what it says, a homomorphism which is bijective. As a consequence two isomorphic groups share many properties, number of elements of a specific order, being abelian, being solvable/nilpotent etc.
Actually isomorphism means two vastly different groups can be understood easily in terms of each other. For example in those days when there were no calculators (for get computers) one used logarithm tables to carry out multiplication. The reason that adding the logarithms and then finding the antilogarithm produces the product of two numbers is because: logarithm gives an isomorphism of the group of positive real numbers under multiplication with the group of all real numbers under addition.
Knowing a newly found  group is isomorphic to a better understood known group  gives us a handle to work with the former.
Many theorems are celebrated because they show two groups occurring in different contexts are isomorphic. (de Rham's theorem in differential geometry/topology, Artin's Reciprocity Law in Algebraic Number Theory)
Easier example: All polynomials with integer coefficients in a single variable is a group under usual addition of polynomials  and is isomorphic to the group of all POSITIVE RATIONAL numbers under usual multiplication of numbers.


TITLE: a trace class operator problem
QUESTION [4 upvotes]: Could someone help me with this
Prove that If $A$ and $B$ are positive trace class operators on a Hilbert space, then so is $A^zB^{(1-z)}$ for a complex number $z$ such that $0 <Re(z)< 1$.
An operator $A$ is called trace class if
$\sum_{n=1}^{\infty}\langle |A|e_n,e_n\rangle < \infty  $
|A| is positive square root of $ A^*A$

REPLY [2 votes]: The trace-class operators are the particular case $p=1$ of the so-called Schatten class $S_p(H)$ which is defined, for $1\leq p<\infty$, by
$$
S_p(H)=\{S\in B(H)\,;\,\|S\|_p=(\mbox{tr}\,|S|^p)^\frac{1}{p}<\infty\}.
$$
that is $|S|^p$ is trace-class.
For $p=2$, these are known as the Hilbert-Schmidt operators. For $1\leq p\leq q$, we have $S_1(H)\subseteq  S_p(H)\subseteq S_q(H) \subseteq K(H)$, where $K(H)$ is the ideal of compact operators. A key fact in this theory is that the noncommutative version of Hölder's inequality holds, including $p$ or $q$ infinite, with $S_\infty(H)=B(H)$. For every $p,q\geq 1$ conjugate exponents, i.e. $\frac{1}{p}+\frac{1}{q}=1$:

$$
\forall S\in S_p(H), T\in S_q(H)\qquad \|ST\|_1\leq \|S\|_p\|T\|_q.
$$
  This says in particular that if $S$ is in $S_p(H)$ and $T$ is in $S_q(H)$, then $ST$ is trace-class. That's the generalization of the well-know case of Hilbert-Schmidt operators $S_2(H)S_2(H)\subseteq S_1(H)$, which is actually an equality, and the ideal property $B(H)S_1(H)B(H)\subseteq S_1(H)$.

Application: let $z=x+iy$ with $0<x<1$, and let $A,B$ be positive trace-class operators, i.e. in $S_1(H)$. Then $|A^z|=A^x$ and $|B^{1-z}|=B^{1-x}$, whence
$$
\left(\|A^z\|_\frac{1}{x}\right)^\frac{1}{x}=\mbox{tr}\left(|A^z|^\frac{1}{x}\right)=\mbox{tr} A=\|A\|_1<\infty$$
and
$$
\left(\|B^{1-z}\|_\frac{1}{1-x}\right)^\frac{1}{1-x}=\mbox{tr}\left(|B^{1-z}|^\frac{1}{1-x}\right)=\mbox{tr} B=\|B\|_1<\infty.
$$
Therefore
$$
\mbox{tr}\big|A^zB^{1-z} \big|=\|A^zB^{1-z}\|_1\leq \|A^z\|_\frac{1}{x}\|B^{1-z}\|_\frac{1}{1-x}=\|A\|_1^x\|B\|_1^{1-x}=(\mbox{tr} A)^x(\mbox{tr} B)^{1-x}<\infty.
$$
In particular, $A^zB^{1-z}$ is trace-class.


TITLE: Differential Equation: Modifying Particular Solution
QUESTION [2 upvotes]: For a particular solution of $$ y^{(4)} - y''' - y'' + y' = t^2 + 4 + te^t$$
The solution to the homogeneous solution is represented by the characteristic equation:
$$ r(r-1)^2 (r+1)$$
so the solutions of the homogeneous solution are:  $ 0,-1$ and $1$ with multiplicity 2.
So the homogeneous solution will have the form:
$$ y_h = Ae^t + Bte^t + Ce^{-t} + D$$
As for the particular solution, I think it will have the form of
$$ y_p = A + Bt + Ct^2 + D + Ee^t + Ete^t$$
and then modify it so that it doesn't overlap with $y_h$.
I am not sure how to go from here. (The question is merely looking for the "form" of the particular solution)

REPLY [3 votes]: $y_p$ can have the form
$$ y_p = A_1 t + A_2 t^2 + A_3 t^3 +  B_1 t^2 e^{t} + B_2 t^3 e^{t}. $$
Here is a general technique tells you how to find $y_p$.
Added: We will apply the general technique to find $y_p$. To annihilate the right hand side of the ode, we apply the operator $D^3(D-1)^2$, where $D=\frac{d}{dt}$, to both sides of the ode which results in
$$ D^4(D-1)^4(D+1) y(t)=0, $$
since $D^3 (t^4+4)=0$ and $(D-1)^2 te^{t} = 0 $.
Now, the above is a new homogenous ode with the the solution,
$$ y_h = c_0 + c_1 t + c_2 t^2 + c_3 t^3 + c_5 e^{t}+ c_6 te^{t}+ c_7 t^2 e^{t} + c_8 t^3e^{t} + c_9 e^{-t} $$
$$ = (c_0 + c_5e^{t}+c_6te^{t} + c_9e^{-t})+ c_1 t + c_2 t^2 + c_3 t^3 + c_7 t^2 e^{t} + c_8 t^3e^{t}.$$
You can see in the above equation what is between the brackets corresponds to $y_h$ of the old ode
$$ y_{h_{o}} = Ae^t + Bte^t + Ce^{-t} + D. $$
So, we can conclude our $y_p$ to be
$$ y_p = A_1 t + A_2 t^2 + A_3 t^3 +  B_1 t^2 e^{t} + B_2 t^3 e^{t}.$$
Note that, the above form will do, but still you can choose the more general form which is based on the above form
$$ y_p = A_0+A_1 t + A_2 t^2 + A_3 t^3 +  (B_0+ B_1  t + B_2 t^2 + B_3 t^3) e^{t} $$


TITLE: If $a_n \leq cb_n$, Then if $\sum b_n$ converges $\sum a_n$ converges.
QUESTION [0 upvotes]: Suppose that $\sum a_n$ and $\sum b_n$ are series of positive terms and that there is a positive number $c$ s.t. $a_n \leq cb_n \forall n$. Show that if $\sum b_n$ converges, then $\sum a_n $ converges.
Is my proof sufficient?
Pf
Assume $\sum a_n$ and $\sum b_n$ are series of positive terms and that there is a positive number $c$ s.t. $a_n \leq cb_n \forall n$
Suppose $\sum cb_n$ converges.
Then $T_n = cb_1 + ... + cb_n = c(b_1 + ... + b_n)$.
We have, $S_n = a_1 + ..... + a_n \leq c(b_1 + ... +b_n) = T_n$
Thus, $S_n$ is bounded monotone sequence.
Therefore $\sum a_n$ converges

REPLY [1 votes]: What you really need to show is that if $\sum b_n$ converges, then $\sum cb_n$ also converges.
Then, use the rest of your proof that compares the partial sums of $a_n$ against the partial sums of $cb_n$ with the inequality.
Now, you need to take a limit of both sides of the inequality to conclude.
At present, you are only saying something about the finite case, which is not enough -- it is the tail end of the sequence that matters, not earlier, possibly "transient behavior" of the sequence.


TITLE: Extrema of two variable function
QUESTION [1 upvotes]: Find extrema of $f(x,y)=x^2-xy+y^2$ from set $M=\{ [x,y] \in \mathbb{R}^2;|x|+|y|\le1\}$

I am solving this kind of problems for the first time and I am not sure what I am doing, what I have got is:

M is compact set(it is bounded and it is closed) and $f$ is continuous, therefore $f$ has maximum and minimum in set M

Firstly I tried to solve it for $Int(M)=\{[x,y] \in \mathbb{R}^2;|x|+|y|<1\}$
$$\nabla f(x,y)=(2x-y,2y-x)=0 \Leftrightarrow (x,y)=(0,0)$$
Which is first "suspicious" point for extreme, but then I got stuck with the the border of my set $B(M)=\{ [x,y] \in \mathbb{R}^2;|x|+|y|=1\}$. I suppose it is possible either without Lagrange multiplier( in somewhat elementary way) or with it, which we just learnt in class recently, and I don't know how to apply it here.

REPLY [1 votes]: The region $ \ |x| + | y | ≤ 1 \ $ is a square with its interior bounded by the four lines $ \ (\pm x) + (\pm y) \ = \  1 \ \ .  $  The function $ \ f(x,y) = x^2 - xy + y^2 \ $ likewise has symmetry about the origin.  We may then expect that two extremal values for the function would occur, one for the points $ \ (\pm x , \pm x) \ $ at which the value of the function is $ \ x^2 - x·x + x^2 \ = \ x^2 \ $ , the other at the points   $ \ (\pm x , \mp x) \ $ with the associated value   $ \ x^2 - x·(-x) + x^2 \ = \ 3x^2 \ \ . $
There are two sets of Lagrange equations here, since there are two pairs of constraint lines:
$ \ \mathbf{x + y = \pm 1  :} \quad 2x - y \ = \ \lambda·1 \ \ , \ \ -x + 2y \ = \ \lambda · 1 \ \ \Rightarrow \ \ \lambda \ = \ 2x - y \ = \ -x + 2y $
$$ \Rightarrow \ \ 3x \ = \ 3y \ \ \Rightarrow \ \ x = y = \pm \ \frac12 $$ $$\Rightarrow \ \ f \left( \pm \ \frac12 \ , \ \pm \ \frac12 \right) \ = \ \left( \pm \ \frac12  \right)^2 - \ \left( \pm \ \frac12  \right)\left( \pm \ \frac12  \right) + \left( \pm \ \frac12  \right)^2 \ = \ \frac14 \ \ ; $$
$ \ \mathbf{x - y = \pm 1  :} \quad 2x - y \ = \ \lambda·1 \ \ , \ \ -x + 2y \ = \ \lambda · -1 \ \ \Rightarrow \ \ \lambda \ = \ 2x - y \ = \ x - 2y $
$$ \Rightarrow \ \ x \ = \ -y \  = \ \pm \ \frac12 $$ $$\Rightarrow \ \ f \left( \pm \ \frac12 \ , \ \mp \ \frac12 \right) \ = \ \left( \pm \ \frac12  \right)^2 - \ \left( \pm \ \frac12  \right)\left( \mp \ \frac12  \right) + \left( \mp \ \frac12  \right)^2 \ = \ \frac34 \ \ . $$
These are two of the values found also by Jason Zimba , which the function attains at the midpoints of the edges of the square.
However, this is a situation that points up the limitations of the Lagrange-multiplier method.  The level curves for the function are rotated ellipses; while the curves $ \ x^2 - xy + y^2 \ = \ \frac14 \ $ and $ \ x^2 - xy + y^2 \ = \ \frac34 \ $ are indeed tangent to opposite sides of the square, the vertices are not covered at all by this approach.  The level curve $ \ x^2 - xy + y^2 \ = \ \mathbf{1} \ $ passes through all four vertices of the square, but cannot be tangent at any of them since the slopes of the sides change discontinuously at the vertices.  The value of the function there needs to be determined in a separate calculation.
Likewise, the "level curve" $ \ x^2 - xy + y^2 \ = \ 0 \ $ is simply the point at the origin, where the value of the function is zero.  So again, the Lagrange method cannot locate this.  We instead find this "critical point" of the function by setting the first partial derivatives $ \ 2x - y \ $ and $ \ -x + 2y \ $ equal to zero; the origin is the intersection of the lines $ \ y = 2x \ $ and  $ \ y = \frac12 x \ $  (as you have shown).  The geometrical situation is represented in the graph below, which may make it clearer why the Lagrange method "misses" the absolute extremal values for this region.


TITLE: fractional chromatic number
QUESTION [1 upvotes]: A fractional chromatic number of a graph $G=(V,E)$, is $min~\sum_{I}y_{I}$, and for every vertex $v\in V$ we have $\sum_{\{I:v \in I\}}y_{I} \ge1$(the condition).
I'm just a little confused on notation.
Let's consider this graph
There are $3$ independent sets: $\{V_{3}\},\{V_{2}\},\{V_{1},V_{4}\}$
So, is a fractional chromatic number of this graph $\le 3$ (sum of all vertices in each independent graph)?
Am I supposed to verify every independent graph if the sum isn't more than $1$, to check the condition ? f.e. for the first set ($\{V_{3}\}$) the sum of all vertices is $1$, and $1 \le 1$, so the condition is met, similarly the next set ($\{V_{4}\}$, and for the last set ($\{V_{1},V_{4}\}$) the sum of all verticies is $\frac{1}{2}+\frac{1}{2}= 1$ and $1 \le 1$ so the condition is met for this and any other set, so indeed, $\chi^{*} \le 3$ ? Is it correct ?

REPLY [0 votes]: There is a variable $y_I$ for each independent set $I$; note that $I$ need not be maximal. Let's call the variables $y_1, y_2, y_3, y_4$, and $y_{14}$. This gives the following optimization problem:
$$\textrm{minimize} \sum_I y_I\\
\begin{align}
\textrm{s.t.}\quad &y_1 + y_{14} &\geq 1\\
& y_2 &\geq 1\\
& y_3 &\geq 1\\
& y_4 + y_{14} &\geq 1\\
\end{align}$$
Clearly $y_2 = y_3 = 1$. For the remaining variables, we'll get the smallest sum when we let $y_{14}$ do all of the work: $y_1=y_4=0$, and $y_{14}=1$. The interpretation of this is that we have no reason to assign different colours to $v_1$ and $v_4$.
So the sum is $3$, and all variables are integer. This corresponds to a $3-$colouring of the graph, which should make sense. (Choose three colours for $\{v_1,v_4\}, v_2$, and $v_3$.) There is no benefit to a fractional colouring for this graph.