Adapting Regularized Low Rank Models for Parallel Architectures

TL;DR

This paper presents a novel smooth non-convex reformulation of regularized low-rank recovery models that leverages parallel architectures, significantly speeding up computations and enabling real-time applications like background subtraction and MRI analysis.

Contribution

It introduces a Burer-Monteiro based reformulation that allows efficient first-order optimization of regularized low-rank models on parallel hardware, with convergence guarantees.

Findings

Achieves an order-of-magnitude speedup over existing RPCA solvers on GPU.

Enables real-time applications such as background subtraction and MRI analysis.

Provides a convergence certificate and suboptimality bounds for the non-convex formulation.

Abstract

We introduce a reformulation of regularized low-rank recovery models to take advantage of GPU, multiple CPU, and hybridized architectures. Low-rank recovery often involves nuclear-norm minimization through iterative thresholding of singular values. These models are slow to fit and difficult to parallelize because of their dependence on computing a singular value decomposition at each iteration. Regularized low-rank recovery models also incorporate non-smooth terms to separate structured components (e.g. sparse outliers) from the low-rank component, making these problems more difficult. Using Burer-Monteiro splitting and marginalization, we develop a smooth, non-convex formulation of regularized low-rank recovery models that can be fit with first-order solvers. We develop a computable certificate of convergence for this non-convex program, and use it to establish bounds on the suboptimality of any point. Using robust principal component analysis (RPCA) as an example, we include numerical experiments showing that this approach is an order-of-magnitude faster than existing RPCA solvers on the GPU. We also show that this acceleration allows new applications for RPCA, including real-time background subtraction and MR image analysis.

Figures

Equations

$$\underset{L,S}{\min }\mathrm{rank}(L)+\text{card}(S),\text{subject to:}\frac{1}{2}\Vert \mathcal{A}(L+S)-b{\Vert }_2^2\le ϵ,$$

$$\underset{L,S}{\min }\Vert L{\Vert }_{\ast }+\Vert S{\Vert }_1,\text{subject to:}\frac{1}{2}\Vert \mathcal{A}(L+S)-b{\Vert }_2^2\le ϵ$$

$$\underset{L,S}{\min }\Vert L{\Vert }_{\ast }+\mathcal{L}(b,\mathcal{A}(L+S))+r(S),$$

$$\underset{L,S}{\min }{\lambda }_L\Vert L{\Vert }_{\ast }+\frac{1}{2}\Vert \mathcal{A}(L+S)-b{\Vert }_2^2+{\lambda }_S\Vert S{\Vert }_1,$$

$$\underset{S}{\min }\left\{{\lambda }_L\Vert L{\Vert }_{\ast }+\frac{1}{2}\Vert \mathcal{A}(L+S)-b{\Vert }_2^2+{\lambda }_S\Vert S{\Vert }_1\right\}={\lambda }_L\Vert L{\Vert }_{\ast }+{\rho }_{{\lambda }_S}(\mathcal{A}(L)-b),$$

$$\phi :L\mapsto \underset{S}{\min }\left\{\mathcal{L}(\mathcal{A}(L+S)-b)+{\lambda }_S\cdot r(S)\right\},$$

$$\underset{L}{\min }{\lambda }_L\Vert L{\Vert }_{\ast }+\phi (L).$$

$$\Vert L{\Vert }_{\ast }\equiv \underset{L=U{V}^T}{\inf }\frac{1}{2}(\Vert U{\Vert }_F^2+\Vert V{\Vert }_F^2).$$

$$\underset{U,V}{\mathrm{argmin}}\phi (U{V}^T)+\frac{{\lambda }_L}{2}(\Vert U{\Vert }_F^2+\Vert V{\Vert }_F^2).$$

$$\underset{U,V}{\min }\mathcal{L}(U{V}^T,X)+r(U)+\tilde{r}(V).$$

$$g(t{L}_1+(1-t)L_2)$$

$$\le f(t{L}_1+(1-t)L_2,t{S}_1+(1-t)S_2)$$

$$\le tf(L_1,S_1)+(1-t)f(L_2,S_2).$$

$$g(t{L}_1+(1-t)L_2)$$

$$=tg(L_1)+(1-t)g(L_2).$$

$$\partial \phi (\cdot )={\mathcal{A}}^{\ast }\nabla \mathcal{L}\left(\mathcal{A}(\cdot )-b-\mathcal{A}(S^{\star })\right),$$

$$\phi (L)=\underset{S}{\min }\left\{\frac{1}{2}\Vert \mathcal{A}(L)+\mathrm{vec}(S)-b{\Vert }_2^2+{\lambda }_S\Vert S{\Vert }_1\right\},$$

$$S^{\star }={\mathrm{prox}}_{{\lambda }_S\Vert \cdot {\Vert }_1}(\mathcal{A}(L)-b).$$

$$\nabla \phi (L)={\mathcal{A}}^{\ast }\left(\mathcal{A}(L)+\mathrm{vec}(S^{\star })-b\right).$$

$$\underset{L}{\min }{\lambda }_L\Vert L{\Vert }_{\ast }+\phi (L)\text{subject to:}\mathrm{rank}(L)\le k.$$

$$\underset{X⪰0}{\min }f(X),\text{such that}\mathrm{rank}(X)\le k,$$

$$\underset{P}{\min }g(P)\stackrel{\text{def}}{=}f(P{P}^T).$$

$$Z=\left[\begin{array}[]{c}U\\ V\\ \end{array}\right]\left[\begin{array}[]{c}U\\ V\\ \end{array}\right]^{T}=\left[\begin{array}[]{cc}UU^{T}&L\\ L^{T}&VV^{T}\\ \end{array}\right].$$

$$f(Z)=\frac{{\lambda }_L}{2}\mathrm{trace}(Z)+\phi (\mathcal{R}(Z)).$$

$$\underset{Z⪰0}{\min }f(Z),\text{such that}\mathrm{rank}(Z)\le k.$$

$$\mathcal{E}\in \partial F(L)$$

$$F(L)-F(L^{\star })\le ⟨\mathcal{E},L-L^{\star }⟩\le \Vert \mathcal{E}{\Vert }_p\cdot \Vert L-L^{\star }{\Vert }_d$$

$$0\in\partial\big{(}\|L\|_{*}+\varphi(L)\big{)}.$$

$$0\in\partial\big{(}\|L\|_{*}+\varphi(L)\big{)}\quad\Leftrightarrow\quad 0\in\big{(}\partial\|L\|_{*}\big{)}+\big{(}\partial\varphi(L)\big{)}.$$

$$L=\widetilde{U}\Sigma\widetilde{V}^{T}=\left[\begin{array}[]{c c}U_{1}&U_{2}\end{array}\right]\left[\begin{array}[]{c c}\Sigma_{1}&0\\ 0&0\end{array}\right]\left[\begin{array}[]{c}V_{1}^{T}\\ V_{2}^{T}\end{array}\right]$$

$$X\in\partial\|L\|_{*}\quad\text{if and only if}\quad\exists W^{\prime}\,\in\mathbb{R}^{m-r\times n-r}\,\text{s.t.}\;\widetilde{U}^{T}X\widetilde{V}=\left[\begin{array}[]{c c}I&0\\ 0&W^{\prime}\end{array}\right],\quad\|W^{\prime}\|\leq 1.$$

$$\Vert \mathcal{E}{\Vert }_F^2$$

Full text

Using Burer-Monteiro splitting and marginalization, we develop a smooth, non-convex formulation of regularized low-rank recovery models that can be fit with first-order solvers. We develop a computable certificate of convergence for this non-convex program, and use it to establish bounds on the suboptimality of any point. Using robust principal component analysis (RPCA) as an example, we include numerical experiments showing that this approach is an order-of-magnitude faster than existing RPCA solvers on the GPU. We also show that this acceleration allows new applications for RPCA, including real-time background subtraction and MR image analysis.

1 Introduction

Low-rank matrix decompositions are an effective tool for large-scale data analytics. Background subtraction, facial recognition, document indexing, and collaborative filtering all use low-rank matrix recovery algorithms (RPCA, ). However, a low-rank structure is only a part of the story — signals often exhibit additional structure. Regularized low-rank recovery models seek to decompose a matrix $X\in\mathbb{R}^{m\times n}$ into the sum of a low-rank component $L\in\mathbb{R}^{m\times n}$ and another structured matrix $S\in\mathbb{R}^{m\times n}$, where $S$ could be sparse, clustered, non-negative, or have some other useful attribute (see e.g. GLRM ).

Searching for sparse and low-rank models gives rise to combinatorial formulations, which are intractable. Fortunately, relaxations often work well, and are amenable to modern optimization techniques. Consider finding an approximate decomposition of a data matrix $X$ into a sum $X=L+S$, where $S$ is sparse and $L$ low-rank. We are given partial linear measurements $b=\mathcal{A}(X)$. We can formulate the problem

[TABLE]

where the cardinality function $\textnormal{card}(\cdot)$ counts the number of nonzero entries of $S$, and the parameter $\epsilon$ accounts for noise or model inconsistencies. Rank is analogous to sparsity, since it is precisely the $\textnormal{card}(\cdot)$ function applied to the singular values of the argument. Problem (1) is NP-hard (NatarajanSparseNP, ) and intractable for large problem sizes. However, both terms admit non-smooth convex relaxations. $\textnormal{card}(\cdot)$ can be replaced with the $\ell_{1}$-norm, and the rank function can be replaced with the nuclear norm, $\|\cdot\|_{*}$, which is equal to the sum of the singular values of the argument:

[TABLE]

When the parameter $\epsilon$ is set to zero, (2) is called robust PCA (RPCA), and can recover the lowest-rank $L$ and the sparsest $S$ under mild conditions. Setting $\epsilon>0$ in (2) is more suitable for most applications, since it allows for noise in the measurements; the problem is then called stable principal component pursuit (SPCP) (RPCA_algo_Wright, ; StablePCP, ). SPCP stably recovers $L$ and $S$ with error bounded by $\epsilon$ (StablePCP, ). Both RPCA and SPCP are convex problems, much easier to solve than the original problem (1). While (2) is an important problem, it is just one example of regularized low-rank models.

Problem class.

We are interested in general regularized low-rank approaches:

[TABLE]

where $\mathcal{L}$ is a loss function and $r$ is a convex regularizer. $\mathcal{L}$, which can be infinite-valued, measures how well $L+S$ agrees with measured values of $X$. In RPCA, $\mathcal{L}(v)=\delta_{b}(v)$ is the indicator function ensuring we match observed data, while in SPCP, $\mathcal{L}(v)=\delta_{2\sqrt{\epsilon}\mathbb{B}}(v-b)$ ensures we are close to $b$. $\mathcal{L}$ can also be finite valued — for example, a modified SPCP is given by solving

[TABLE]

where $\lambda_{L},\lambda_{S}$ are tuning parameters, and $\mathcal{L}(\cdot)=\frac{1}{2}\|\cdot\|_{2}^{2}$.

The regularizer $r(S)$ is a convex penalty that promotes the desired structure of $S$. The $\ell_{1}$ norm $r(S)=\|S\|_{1}$ promotes sparsity in $S$; other important examples include the ordered-weighted $\ell_{1}$ norm (OWL norm) SLOPE ; davis2015 ; OWLnorm , and the elastic net (ElasticNet, ), which both enforce sparsely correlated observations (clusters) in $S$. Other convex regularizers promote known group-sparse and hierarchical structures in $S$; see BachPenalties for a survey.

Marginalization in $S$.

In all of the applications we consider, $r(S)$ is simple enough that we can efficiently minimize (3) in $S$ for a given $L$. In the motivating SPCP example (4), we have

[TABLE]

where $\rho_{\lambda_{S}}$ is the Huber function (for a detailed derivation, see e.g. LevelSet ). This motivates the definition

[TABLE]

where $\varphi$ is smooth (see Section 3.2) with value and derivative either explicitly available as in the case of SPCP, or efficiently computable. After marginalizing out $S$, (3) becomes

[TABLE]

Proximal gradient methods are applicable, but require computing the expensive and communication-heavy singular value decomposition (SVD) at each iteration. To ameliorate this bottleneck, it is possible to parallelize the SVD step using communication-avoiding algorithms. We develop such an approach in Section 6.1 to use as a performance baseline.

However, the focus of this paper is to use a non-convex, factorized representation of the nuclear norm that avoids the SVD entirely. This reformulated program offers more speedup through parallelization than communication-avoiding approaches, and the introduced non-convexity can be easily controlled. The smooth structure of $\varphi$ allows a computable optimality certificate for the non-convex case, which is developed and explored in Section 3.

Non-convex reformulation.

To avoid the SVD entirely, we consider non-convex reformulations for regularized low-rank recovery models, using the factorization $L=UV^{T}$, where $U\in\mathbb{R}^{m\times k},\ V\in\mathbb{R}^{n\times k}$, and $k\ll m,n$. This factorization imposes the rank-constraint $\mathrm{rank}(L)\leq k$. In addition, we can further penalize rank using the representation (see e.g. BachSparseFactor ; Recht07 ; MMMF2005 )

[TABLE]

This equivalence only holds if the ranks of $U$ and $V$ are large enough, and makes it possible to maintain a low-rank $L=UV^{T}$ without requiring SVDs. The non-convex analogue to (6) is given by

[TABLE]

Roadmap.

We provide a survey of related work in Section 2. We consider two theoretical issues related to the smooth reformulation (8) in Section 3. First, the factorization $L=UV^{T}$ gives a non-convex program that could introduce suboptimal stationary points and local minima. Second, it is important to understand the properties of the map $\varphi$ obtained by marginalizing over $S$ in the context of specific optimization algorithms. In Section 4, we study heuristic initialization methods, which are important for non-convex programs. We present detailed numerical studies, comparing the non-convex method with competing alternatives, in Section 5. In this section, we also demonstrate the effectiveness of our non-convex formulation of low-rank recovery models for real-time video processing. There are several approaches to nuclear norm minimization that focus on accelerating the SVD step, either by reducing its communication costs or by replacing it with a simpler subproblem, as in the Frank-Wolfe method. In Section 6, we develop two such methods for general regularized low-rank recovery models, and we compare these to our non-convex approach. Our numerical experiments illustrate that our non-convex program is able to solve regularized low-rank recovery problems faster than existing solvers — in many cases, by an order of magnitude — without sacrificing accuracy in the solution. We end with final thoughts in Section 7.

2 Prior Work and Contributions

Several authors have used alternating minimization to solve non-convex, factorized matrix completion models (GAGG13, ; HardtAltMin, ; MatCompAltMin, ; KeshPhD, ; RegForMatComp, ; NoisyMatComp, ). This line of research is supported by the existence of conditions that guarantee fast convergence to a globally optimal solution (GAGG13, ; HardtAltMin, ; RegForMatComp, ; NoisyMatComp, ), sometimes at a linear rate (MatCompAltMin, ; GLRM, ). These models are a special class of (8), when no regularizer is present (i.e., $r(S)\equiv 0$).

In GLRM , the authors review the use of alternating minimization to solve regularized PCA problems of the following form:

[TABLE]

Matrix completion models and RPCA can adopt the split-form in (9), but SPCP and other regularized models cannot. Our analysis applies to more general regularizers and offers an approach using first-order solvers rather than alternating minimization.

The authors of Shen2012 develop a split RPCA program and solve it using alternating minimization as well. While their technique can be an order of magnitude faster than competitive convex solvers (such as the inexact augmented Lagrangian method (IALM) from IALM ), its performance suffers when the magnitude of the sparse component is large compared to the magnitude of the low-rank component.

In Procrustes , the authors consider provable recovery of a low-rank solution to a system of linear equations using a non-convex, factorized program. In a similar vein, the work of CaramanisRPCA extends these ideas to investigate a non-convex formulation of robust PCA that can be solved with gradient descent.

Another non-convex approach to RPCA is considered in NonConRPCA . In contrast to the research discussed previously, their method does not use a low-rank factorization of $L$. Instead, these authors develop an algorithm of alternating projections, where each iteration sequentially projects $L$ onto a set of low-rank matrices and thresholds the entries of $S$. Their technique also comes with recovery and performance guarantees.

We develop a general non-convex formulation of regularized low-rank recovery models that can be solved using first-order optimization algorithms. We address problems associated with local minima and spurious stationary points that accompany non-convexity, by showing that all local minima of our model have the same objective value, and by providing a certificate to prove that a given solution is not a spurious stationary point. We also show that our method is particularly well-suited for the GPU.

3 Theoretical Considerations

We discuss three issues related to the approach (8). First, we study the convexity and differentiability of the marginal function (6) in Sections 3.1 and 3.2. Next, moving from (6) to (8) by factoring $L$, we transform a convex problem into a non-convex one. This creates two issues: an implicit rank constraint $\operatorname{rank}(L)\leq k$ that is not present in (8), and potential for local minima and spurious stationary points in the non-convex variant that do not correspond to a solution of the convex problem. We address these issues in Section 3.3.

Section 3.4 introduces a computable certificate that can be used to check whether the non-convex model has converged to a global solution. Using the smoothness of $\varphi$, we derive a simple optimality certificate that can be tracked while optimizing (8) to identify either when the implicit rank of the factors is too small or when we have converged to a minimizer of (6).

3.1 Convexity of $\varphi$

The convexity of the marginal function (6) is a well-known result in convex analysis, see e.g. (VariationalAnalysis, , Proposition 2.22). We include a short proof of this fact.

Lemma 1**.**

Let $f:\mathbb{R}^{m\times n}\times\mathbb{R}^{m\times n}\to\mathbb{R}$ be a convex function, and define $g(L)=\min_{S\in\mathbb{R}^{m\times n}}f(L,S)$ for $L\in\mathbb{R}^{m\times n}$. The function $g$ is convex.

Proof.

For all $L_{1},L_{2},S_{1},S_{2}\in\mathbb{R}^{m\times n}$ and $t\in(0,1)$,

[TABLE]

As this inequality holds for arbitrary $S_{1}$ and $S_{2}$, it holds for any minimizing $S^{\star}$:

[TABLE]

∎

It follows immediately that $\varphi$ as defined in (5) is convex when $\mathcal{L}$ and $r$ are convex.

3.2 Smoothness of $\varphi$

When $\mathcal{L}$ is the least-squares error and $\mathcal{A}$ is an injective linear operator, Lipschitz differentiability of $\varphi$ is immediate from the strong convexity of $\mathcal{L}$. Support for this claim is given, for example, in CombettesBook , where the relevant proposition for the specific case where $\mathcal{A}$ is a reshaping operator is as follows:

Proposition

(CombettesBook**, **, Prop. 12.29)****: Let $r$ be a proper lower-semicontinuous convex function, and let ${}^{\mu}r(\mathcal{A}(L)-\mathrm{vec}(S)-b)=\inf_{S}\ r(S)+\frac{1}{2\mu}\|\mathcal{A}(L)-\mathrm{vec}(S)-b\|_{2}^{2}$ be the Moreau envelope of $f$ with parameter $\mu$. Then ${}^{\mu}r$ is differentiable and its gradient is $\mu^{-1}$-Lipschitz continuous.

More generally, $\varphi(L)$ is differentiable as long as $r(S)$ has a unique minimizer and the loss function $\mathcal{L}$ is convex and smooth. By (VariationalAnalysis, , Theorem 10.58), we have

[TABLE]

where $S^{\star}$ is any minimizer of (4). The subdifferential is a singleton if and only if the minimizing $S^{\star}$ is unique, and a unique minimum is guaranteed if the sum $\mathcal{L}(\mathcal{A}(L+S)-b)+r(S)$ in (4) is strictly convex.

Using SPCP as an example, we have

[TABLE]

where sparse outliers $S$ are fully observed. The objective is strongly convex in $S$, and has the unique minimizer

[TABLE]

Therefore, $\varphi(L)$ is smooth, with

[TABLE]

We used (10) rather than the parametric form of $\varphi(L)$ to obtain its value and gradient. Using (7) to replace the nuclear norm and $\varphi$ to eliminate the non-smooth regularizer, we now have a smooth optimization problem amenable to first-order methods.

For minimizing smooth objectives, quasi-Newton methods are often faster than algorithms that require only first-order smoothness, such as gradient-descent. The convergence rates of quasi-Newton methods depend on second-order smoothness, which does not hold in problems such as SPCP, since e.g. the Huber function is only $\mathcal{C}^{1}$. However, empirical results (including those presented in Section 5) suggest that these methods are effective.

3.3 Rank and Local Minima

The factorized problem (8) is not equivalent to (6) because of the implicit rank constraint. However, we can show that (8) is equivalent to the rank-constrained problem

[TABLE]

This follows immediately from (7). In order for the split problem to recover the solution to (1), it must be initialized with a $k$ larger than the rank of the minimizing $L$. However, a larger $k$ slows computation, so it must be chosen with care. This issue is considered in more depth in Section 4.

We want to be sure that any local minimum of (8) corresponds to a (global) minimum of (6). The following theorem combines ideas from BurerMonteiro2005 and Aravkin2014 to show this holds provided that $k$ is larger than the rank of the minimizer for (6).

Theorem 1**.**

Consider an optimization problem of the following form:

[TABLE]

where $X\in\mathbb{R}^{n\times n}$ is a positive semidefinite real matrix, and $f$ is a lower semi-continuous function mapping to $[-\infty,\infty]$ and has a non-empty domain over the set of positive semi-definite matrices. Using the change of variable $X=PP^{T}$, take $P\in\mathbb{R}^{n\times k}$, and consider the problem

[TABLE]

Let $\bar{X}=\bar{P}\bar{P}^{T}$, where $\bar{X}$ is feasible for (12). Then $\bar{X}$ is a local minimizer of (12) if and only if $\bar{P}$ is a local minimizer of (13).

The proof of Theorem 1 is deferred to Appendix A. Using the SDP formulation of the nuclear norm presented in fazel2001rank , our problem can be recast as a semi-definite program so that we can apply Theorem 1. Define

[TABLE]

The matrix $Z$ is positive semi-definite, and has form $Z=PP^{T}$, with P=\left[\begin{array}[]{c}U\\ V\\ \end{array}\right]. Let $\mathcal{R}(\cdot)$ be the function that extracts the upper-right block of a matrix (so that $\mathcal{R}(Z)=L$), and let

[TABLE]

Using equation (3), we now see that the rank-constrained problem is equivalent to

[TABLE]

Applying Theorem 1, we can now be assured that $\overline{P}$ is a local minimizer to the split program if and only if $\overline{Z}$ is local minimizer of the original problem. This is equivalent to the statement that the point $(\overline{U},\overline{V})$ is a local minimizer of the split problem if and only if $\overline{L}=\overline{U}\ \overline{V}^{T}$ is a local minimizer of the original, rank-constrained program.

3.4 Certificate of Convergence

Although Theorem 1 asserts that the split problem and the rank-constrained problem have the same local minima, the rank-constrained problem is itself non-convex, so we have not guaranteed that every stationary point of the split problem solves the convex program. Using an approach that builds on GLRM , we develop a certificate to check whether a recovered solution to (6) corresponds to a spurious stationary point when $\varphi$ is convex (in particular when $\mathcal{L}$ and $r$ are convex). This technique can also be used on the convex formulation as a method to check the distance to optimality.

We base our certificate on the following analysis. Let $F$ be any proper convex function, then Fermat’s rule states that the set of minimizers of $F$ are points $L$ such that $0\in\partial F(L)$, where $\partial(\cdot)$ denotes the subdifferential. One could take a candidate point $L$, construct $\partial F(L)$, and if it contains [math], conclude that $L$ is the minimizer; this is not realistic for large-scale problems. Instead we find

[TABLE]

where $\mathcal{E}$ is small. Let $L^{\star}$ denote any minimizer of $F$, then from the definition of a subgradient, $F(L^{\star})\geq F(L)+\langle\mathcal{E},L^{\star}-L\rangle$. Then

[TABLE]

for any pair of primal and dual norms $\|\cdot\|_{p}$ and $\|\cdot\|_{d}$. We discuss bounding $\|\mathcal{E}\|_{p}$ and $\|L-L^{\star}\|_{d}$ in the next two subsections.

Finding an approximate zero $\mathcal{E}$

The key to the certificate is analyzing the marginalized problem (6) in terms of $L$, rather than the problem (3) in terms of $(L,S)$. The point $L=UV^{T}$ is the optimal solution of (6) if and only if

[TABLE]

We have set $\lambda_{L}=1$ for convenience, or alternatively one can absorb $\lambda_{L}^{-1}$ into $\varphi$. Since both the nuclear norm and $\varphi$ are proper lower semi-continuous and the intersection of the interior of their domains is non-empty (in particular they are finite valued), then by (CombettesBook, , Cor. 16.38), we have that

[TABLE]

Both of these subdifferentials are computable. Let

[TABLE]

be the (full) SVD of $L$, where $U_{1}\in\mathbb{R}^{m\times r}$ and $V_{1}\in\mathbb{R}^{n\times r}$. The subdifferentials of the nuclear norm at $L$ comprises matrices of the form $X=U_{1}V_{1}^{T}+W$, where $U_{1},V_{1}$ contain the left and right singular vectors of $L$ that correspond to non-zero singular values,111The orthogonal matrices $U_{1}$ and $V_{1}$ should not be confused with the variables $U$ and $V$ that are used as a factorization of $L$, as $U$ and $V$ are not necessarily orthogonal. In fact, $U_{1}$ and $V_{1}$ can be efficiently computed from the factorization $L=UV^{T}$ by taking the QR-decompositions $U=Q_{U}R_{U}$ and $V=Q_{V}R_{V}$ so $L=Q_{U}(R_{U}R_{V}^{T})Q_{V}^{T}$. Perform a SVD on the small inner matrix to write $(R_{U}R_{V}^{T})=U_{R}\Sigma V_{R}^{T}$, and hence $\Sigma$ are the nonzero singular values of $L$ and $U_{1}=Q_{U}U_{R}$ and $V_{1}=Q_{V}V_{R}$ are the corresponding left and right singular vectors.

and $W$ satisfies the conditions $U_{1}^{T}W=0$, $WV_{1}=0$ and $\|W\|_{2}\leq 1$. Equivalently,

[TABLE]

We can find an explicit matrix $D\in\partial\varphi(L)$ using (10), and it is unique under the smoothness conditions described in Section 3.2.

Altogether, we can guarantee a point $\mathcal{E}\in\partial\|L\|_{*}+\partial\varphi(L)$ such that

[TABLE]

In the first equation, the minimum is achieved since the squared-norm is continuous and the set is compact, and (21) is used in the third equation. Each term of the fourth equation can be calculated explicitly, with the last term obtained via projection onto the unit spectral-norm ball, which requires the SVD of the small matrix $U_{2}^{T}DV_{2}$. This bound on $\|\mathcal{E}\|_{F}$ can then be used in (18).

Most terms above can be computed efficiently, in the sense that $U_{2}$ and $V_{2}$ never need to be explicitly computed (which is most important when $r$ is small, since then $U_{2}$ and $V_{2}$ have $m-r$ and $n-r$ columns, respectively), and therefore the computation consists only of matrix products, thin-QR factorizations, and an $r\times r$ SVD, leading to a complexity of $\mathcal{O}(rmn+r^{2}\cdot(m+n)+r^{3})$. For example, we compute

[TABLE]

and $\|U_{2}^{T}DV_{1}\|_{F}^{2}$ is computed analogously. The final term requires an unavoidable $(m-r)\times(n-r)$ SVD factorization, so for large $m,n$ the certificate may be computed as a final check rather than at each iterate.

Bounding the distance to the feasible set

We seek a bound on $\|L-L^{\star}\|_{d}$ in an appropriate norm, which together with (18) gives us a bound on the objective function difference $F(L)-F(L^{\star})$. Since the bound on $\mathcal{E}$ is in the Frobenius norm, ideally we set $\|\cdot\|_{d}=\|\cdot\|_{F}$, but bounds in any norm will work using $\|L\|_{F}\leq\sqrt{\text{rank}(L)}\|L\|$ and $\|L\|_{F}\leq\|L\|_{*}$.

Letting $F(L)=\|L\|_{*}+\varphi(L)$, we first bound $F(L^{\star})$ by computing $F(L)$ for explicit choices of $L$. To be concrete, in this section we assume $r(S)=\|S\|_{1}$ with $\lambda_{S}=1$, and $\mathcal{L}(\cdot)=\frac{1}{2}\|\cdot\|^{2}$. Choosing $L$ such that $\mathcal{A}(L)=b$ means that $S=0$ in the definition of $\varphi$ (5), and hence $F(L)=\|L\|_{*}$. Choosing $L=0$ and $S=0$ gives the bound $F(L)\leq\frac{1}{2}\|b\|^{2}$. A third choice is to explicitly compute $F(L)$ at all iterates in the algorithm and record the best.

Denoting $F^{\star}=\min_{L}F(L)$ and using $F^{\star}\leq F_{\text{bound}}$, non-negativity of $\varphi$ immediately implies $\|L^{\star}\|_{*}\leq F_{\text{bound}}$. Hence $\|L-L^{\star}\|_{F}\leq\|L\|_{F}+\|L^{\star}\|_{F}\leq\|L\|_{F}+\|L^{\star}\|_{*}\leq\|L\|_{F}+F_{\text{bound}}$, and $\|L\|_{F}$ is explicitly computable.

Results

Figure 1 shows how the distance to the optimal subgradient decreases over time for the convex solver LagQN, the non-convex Split-SPCP with a rank bound that is too strict, and Split-SPCP with a rank bound large enough to reach the global optimum. When the rank bound is too restrictive, Split-SPCP cannot recover the optimal $L$. Measuring the distance to the optimal $L$ reveals this. Figure 1 shows that the distance plateaus far from zero when the rank of $L$ is bounded above by 30, and the rank of the optimizing $L$ is 58. In practice, this measure can be used to indicate that the rank bound should be increased.

4 Initialization

The non-convex formulation (8) is sensitive to the initial point, and requires a bound on the rank of $L$. In this section, we present heuristics for choosing a good initial point and determining an appropriate rank-bound.

4.1 Dynamically Increasing $k$

Figure 2 demonstrates the sensitivity of the Split-SPCP program to the factor rank $k$. These tests were performed on the surveillance video data described in Section 5, where the true rank of the low-rank component is 58. We see that when $k\leq\operatorname{rank}(L)$, Split-SPCP does not converge to the correct solution, but if $k$ is much larger than the rank of the minimizing $L$, then the computation is slowed significantly.

Because our solver is oblivious to the individual dimensions of $U$ and $V$, columns can be added to both matrices on the fly. Dynamically updating the rank bound can help when (8) is initialized with a $k$ that is too small. Suppose a solver reaches some convergence criterion at iteration $i$. To see if this is a solution to (6), we add a single column to $U$ and $V$:

[TABLE]

and observe whether this rank-one allows for a lower objective value or certificate value. Dynamically increasing $k$, even aggressively, is more efficient than overestimating $k$ from the start.

4.2 Choosing $U_{0}$ and $V_{0}$

Choosing an initialization point $(U_{0},V_{0})$ has received considerable attention in factorized low-rank recovery problems. Most of this work has centered around matrix completion using alternating minimization, see, for example, HardtAltMin ; MatCompAltMin ; RegForMatComp ; NoisyMatComp ; GLRM . Provably good choices for $U_{0}$ and $V_{0}$ are $U_{0}=U\Sigma^{\frac{1}{2}}$ and $V_{0}=V\Sigma^{\frac{1}{2}}$, where $U$ and $V$ come from the SVD of $\mathcal{A}^{*}(X)$, where $\mathcal{A}$ is a sampling operator and $\mathcal{A}^{*}$ is its adjoint. HardtAltMin showed that these initial points lie in the basin of attraction of globally optimal minima, yielding fast convergence to an optimal solution.

As mentioned in GLRM , it is sometimes not necessary to perform a full SVD to form $U_{0},V_{0}$. Once $k$ is chosen, only a partial SVD is necessary to calculate the first $k$ singular values and vectors, which can be done efficiently using the randomized SVD (rSVD) (structureRandomness, ). Although using the rSVD is significantly faster than a full SVD, the values and vectors it returns are not always accurate, especially when singular values of the underlying data do not decay rapidly.

Figure 3 shows the performance of Split-SPCP using various initial points. The tests were performed on a $5,\!000\times 5,\!000$ matrix $X=L+S$, where $\operatorname{rank}(L)=1000$ and $\textnormal{sparsity}(S)$ = 9.43%. The rank bound was $k=1,\!050$. All $5,\!000$ singular value and vector triples of $X$ were calculated for the SVD test, and $U_{0}$ and $V_{0}$ were formed from the first $k=1,\!050$ triples. For the rSVD test, only the first $1,\!050$ triples were approximated, $U_{0}$ and $V_{0}$ were formed from these.

Figure 3 shows that initializing $U_{0}$ and $V_{0}$ with the first 1,050 singular value and vector triples returned by the full SVD yields the smallest initial error, but the extra computational cost is not worth the better start. The rSVD allows quicker convergence to a similar solution. The two random initializations converge to a stationary point that is not globally optimal.

5 Numerical Experiments

In this section, we present numerical experiments that illustrate the speed of Split-SPCP. Recall that the Split-SPCP program is given by

[TABLE]

where

[TABLE]

Since the objective of Split-SPCP is smooth, any first-order method can be used to solve it. Also, because the communication-heavy SVD step is no longer a limitation, we are motivated to choose a solver and an implementation that are most suited for the GPU.

5.1 Implementation Details

To find a solver that fits our Split-SPCP program well, we compare the performance of several first-order methods as they solve a logistic regression problem on the GPU. Each of the solvers we use in this test are modified from SchmidtCode , using MATLAB’s Parallel Computing Toolbox, which in turn uses MAGMA libraries (Magma, ) to run on the GPU. For all GPU computation, we use a Tesla K40c GPU card with thread block sizes $[1024,1024,64]$, grid sizes $[2.15\times 10^{9},6.6\times 10^{4},6.6\times 10^{4}]$, a clock rate of 745 MHz, and 12.0 GB of total memory. For all of the tests run on the CPU in this section, we compute on an Intel Xeon CPU E5-2630 version 3, using a 64-bit architecture and a clock rate of 2.40GHz. Our version of MATLAB is R2016b.

The time it takes for various solvers to minimize a logistic loss with a Tikhonov regularization term for $10^{4}$ labels is shown in Table 1 below. We measure performance in “computation time”, the time it takes for each solver to come within $10^{-8}$ of the optimal objective value.

We see from our analysis that the L-BFGS solver, storing a small number of iterations in memory, yields the best performance. Of course, the relative performance of different solvers on the GPU depends on the stopping tolerance as well as the specific problem, but there is further evidence in the literature supporting L-BFGS as a competitively efficient solver (WrightBook, ).

5.2 SPCP for Background Subtraction

In ThierryVideoReview and Sobral , the authors provide a review of the fastest algorithms in the family of robust PCA and principal component pursuit for background subtraction. We choose some of the fastest of these to compare to our Split-SPCP algorithm. Each of these methods runs faster on the GPU than the CPU, so we run all of them on the GPU for our comparisons (using MATLAB’s Parallel Computing Toolbox). We find that our non-convex solver is particularly well-suited for the GPU, and outperforms all other solvers in almost every case.

For our background subtraction tests, we use the escalator surveillance video provided by LiData . We want to identify the people in the video while ignoring the moving escalators and the stationary background. This problem is particularly difficult for many low-rank recovery algorithms, such as PCA, because the motion of the escalator is a confounder. SPCP is less sensitive to outliers, so it can overcome this challenge.

To find a reference solution, we hand-tune the parameters $\lambda_{L}$ and $\lambda_{S}$ in the (convex) quasi-Newton Lagrangian SPCP algorithm (LagQN) until we find a qualitatively accurate decomposition. This process produces a reference low-rank component $L_{ref}$ and sparse component $S_{ref}$. The $L_{ref}$ we choose has rank 58, and $S_{ref}$ is 58.03% sparse. The optimal parameters for both Split-SPCP and LagQN are $\lambda_{L}=115$ and $\lambda_{S}=0.825$. We tune the parameters in the other solvers to recover this solution as closely as possible. We measure error as the normalized difference between the objective and the objective at the reference solution. Since the solvers minimize different objectives, we calculate the Split-SPCP objective value at each iteration for every algorithm, and we do not include these calculations in our time measurements. To initialize Split-SPCP, we use the first 100 singular values and vectors from the randomized SVD of the data matrix.

Several of the algorithms in Figure 5 do not converge to the same solution, despite considerable effort in parameter tuning. These algorithms were designed to quickly find approximate solutions to the SPCP problem, and might suffer from the large amount of noise present in the data. The approximate solutions recovered by these algorithms are qualitatively different from the Split-SPCP solution, see Figure 6. We also find that a lower objective value generally corresponds to a qualitatively superior solution.

To quantitatively measure the quality of the solutions, we use the Corrected Akaike Information Criterion, or AIC${}_{\textnormal{\large c}}$. The AIC${}_{\textnormal{\large c}}$ measures the fit of a statistical model to a certain data set. Given a model with $p$ parameters and a log-likelihood function $\ell$, the value of the AIC${}_{\textnormal{\large c}}$ is

[TABLE]

where $m\cdot n$ is the size of the data set and $\ell_{\max}$ is the maximum of $\ell$. The preferred statistical model is the one that minimizes the AIC${}_{\textnormal{\large c}}$. The AIC${}_{\textnormal{\large c}}$ favors models that maximize the likelihood function and penalizes complex models with many parameters, guarding against overfitting. Using AIC${}_{\textnormal{\large c}}$ as a measure of quality avoids both overemphasis of objective values and the inherent ambiguity of visual comparisons.

To compute the AIC${}_{\textnormal{\large c}}$ value, we must formulate SPCP as a statistical model. It is well-known that the least-squares loss term assumes that the data is of the form $X=L+S+Z$, where the entries of $Z$ are i.i.d. Gaussian random variables with $\mu=0$ and variance estimated using the sample variance, $\hat{\sigma}^{2}=\frac{\|X\|^{2}_{F}}{m\cdot n}$. Similarly, the $\ell_{1}$-regularizer assumes that the entries of $S$ are drawn iid from the Laplace distribution with mean 0 and variance estimated as $\hat{b}=\frac{\|S\|_{1}}{m\cdot n}$. The nuclear norm of $L$ is the $\ell_{1}$-norm of its singular values, so its corresponding prior assumes that the singular values of $L$ follow the Laplace distribution. The log-likelihood function is then

[TABLE]

where $\sigma^{2},b$ and $b_{*}$ are computed according to their respective estimator. We must also define the number of parameters of this model, which is equal to its degrees of freedom. Each term in SPCP provides the following degrees of freedom:

[TABLE]

Note that $nnz(S)$ counts the number of non-zero entries in $S$. Finding the degrees of freedom in a rank-$k$ matrix or a sparse matrix are both standard calculations. For the loss term, we use the residual effective degrees of freedom as an estimate of the flexibility that it introduces. This is detailed further in LoaderRegression . The number of parameters $p$ is equal to the total degrees of freedom.

The AIC${}_{\textnormal{\large c}}$ values for the solutions shown in Figure 6 are listed in Table 2. The AIC${}_{\textnormal{\large c}}$ value for the reference solution is listed under “oracle.” Recall that the reference solution was found using the LagQN solver with a tight tolerance. All the other values in the table correspond to solutions that meet comparable tolerances. Also, since the IALM model does not use an $\ell_{1}$-norm regularizer, many of the values in the returned $S$ matrix were small but not exactly zero. To make accurate comparisons, we applied the shrinkage operator to the $S$ matrix returned by IALM to set small values equal to zero before calculating the degrees of freedom. The values for IALM and LMaFit are high because these solvers are not generally robust to large amounts of noise.

With the AIC${}_{\textnormal{\large c}}$ metric as well, we see that Split-SPCP discovers the best solution, and in a much shorter time compared to the other algorithms.

Split-SPCP has an advantage because it parallelizes well on GPU architectures. This is due to the relative simplicity of the L-BFGS algorithm, which avoids complex linear algebraic decompositions that require heavy communication, such as the QR and SVD (see also ParallelLBFGS ). Figure 7 shows that although Split-SPCP on the CPU is slower than LagQN on the CPU, Split-SPCP enjoys enormous speedup when it is implemented on the GPU.

5.3 Synthetic Data

Split-SPCP demonstrates similar performance on synthetic data. For our synthetic-data test, we create two $1,\!000\times 150$ random matrices with each entry drawn from independently from the univariate normal distribution, and define our low-rank reference matrix, $L_{ref}$, as the product of these two matrices. To form the sparse reference matrix, $S_{ref}$, we fill 50% of a $1,\!000\times 1,\!000$ matrix with numbers drawn independently from the univariate normal distribution, and the other entries we set equal to zero. $X$ is then the sum $L_{ref}+S_{ref}$ with added white Gaussian noise, so that $\frac{\|L_{ref}+S_{ref}-X\|_{F}}{\|X\|_{F}}=8.12\times 10^{-5}$. Matrices with these characteristics are often encountered in video processing applications.

We initialize Split-SPCP with the first 200 singular value and vector triples of $X$, so that $k$ is larger than the rank of $L$. These triples are found using the rSVD. The tuning parameters are set to $\lambda_{L}=2.95$ and $\lambda_{S}=0.1$. As before, we measure performance based on the normalized difference from the true Split-SPCP objective value. The results are shown in Figure 8.

We see that FPCP and GoDec discover approximate solutions quickly, but the approximations are not within $10^{-2}$ of the true objective. LMaFit exhibits poor performance. This could be due to the fact that $S$ is about 50% sparse, and LMaFit struggles with problems that have a large sparse component (LMaFit, ). Split-SPCP significantly outperforms all solvers in speed, and outperforms all solvers except for LagQN in accuracy. LagQN and Split-SPCP converge to solutions with the same objective value, albeit Split-SPCP converges much more quickly.

5.4 Dimensional Scaling and Real-Time Video Processing

Background subtraction in video is a quintessential application for low-rank recovery algorithms, but these algorithms are generally unable to perform real-time background subtraction due to their slow implementation. The experiments in Section 5.2 show that Split-SPCP is faster and more accurate than the methods reviewed in ThierryVideoReview and Sobral for background subtraction, and in this section, we use Split-SPCP for the real-time processing of high-resolution videos with standard frame rates. We also explore how Split-SPCP scales as the resolution and frame rate increases.

To measure the performance of Split-SPCP on more diverse video datasets, we use videos from the Background Model Challenge (BMC, ). The experiments shown in Figure 9 consider a video with $240\times 320$ pixels per frame and 32,964 frames, which is substantially larger than the low-resolution escalator video of Section 5.2 with $130\times 160$ pixels per frame and only 200 frames. We unfold the three-dimensional video data into a matrix with each column containing one frame, and we partition the matrix columnwise into blocks comprising a fixed number of frames. With the rank bound fixed at 60, we record the computation time required to perform 200 iterations of Split-SPCP on one block as a function of the block size. Larger block sizes yield better solutions, but they also require more time for SPCP to converge. We find that after 200 iterations, Split-SPCP converges to within $10^{-3}$ of the optimal objective value for all block sizes, and the solution is qualitatively similar to the optimal solution.

Figure 9 shows how Split-SPCP scales as the video resolution and the number of frames per block increase. These two quantities correspond to the number of rows and the number of columns of each block. We increase the resolution by rescaling each frame so that the aspect ratio stays close to constant, and we use interpolation to impute the missing pixel values. It is easy to discern the $\mathcal{O}(n)$ scaling as the number of frames (columns) is increased. Increasing the resolution affects the computation time linearly, but beyond a “critical resolution” of $332\times 446$ pixels per frame, the slope of the linear scaling undergoes a shift in regime. This linear scaling is expected because the bottleneck of Split-SPCP is the matrix-matrix multiply $UV^{T}$, which requires $\mathcal{O}(mnk)$ operations, so with $k$ fixed, the computation time of Split-SPCP grows linearly with the number of data points. The results of the background subtraction are shown in Figure 10.

These experiments suggest conditions under which Split-SPCP can decompose a video in real-time. Assuming a resolution of $240\times 320$ pixels per frame and a frame rate of 24 frames-per-second, which is standard in many applications of video processing, the video must be partitioned blocks of about 600 frames or more. With this partitioning, the algorithm will have finished decomposing the given block before the next block is recorded. This number increases linearly with resolution, following the trend line given in the rightmost plot of Figure 9. Following these guidelines, Split-SPCP can decompose the videos from BMC in real-time without sacrificing quality. An example is shown in Figure 10.

6 A Comparison with Accelerated-SVD Algorithms

There are several algorithms for low-rank matrix recovery that focus on accelerating the SVD step. In this section, we explore two such approaches: a communication-avoiding rSVD and Frank-Wolfe with marginalization. We compare both to our SVD-free approach. Both of these methods preserve convexity, but numerical experiments suggest that our non-convex formulation is superior.

6.1 Communication-Avoiding rSVD

Many problems in data analysis involve matrices that are “tall-skinny”. For example, the surveillance video we consider in Section 5 has dimensions $20,\!800\times 200$, and the psychometric data shown in Figure 12 produces a matrix with 230-times more rows than columns. The prevalence of this tall-skinny structure is due in part to the fact that multi-dimensional data sets must be “unwrapped” into a two-dimensional matrix before matrix-based algorithms can analyze the data, and the canonical unwrapping of a square data set with more than two dimensions leads to a highly rectangular matrix. Tall-skinny matrices pose a challenge for matrix decompositions, especially on computational architectures where communication is costly (such as a GPU) (TSQR, ).

The rSVD is often used to accelerate low-rank recovery models by projecting a data matrix of size $m\times n$ onto a low-dimensional subspace, and performing computations with the resulting smaller matrix of size $m\times k$, with $k\ll n$.222Further details of the rSVD algorithms can be found in halko2011finding The rSVD offers improved performance on the CPU, but it exacerbates the tall-skinny problem, making parallelization difficult even when the input matrix is square.

The effect of this increased communication is notable: on the CPU, the rSVD step in LagQN accounts for about 55% of the algorithm’s total running time, but when this solver is run on the GPU, the rSVD accounts for about 90% of SPCP’s total running time, as shown in Figure 11. The TSQR algorithm of TSQR can provide a convex option for accelerating low-rank recovery models. While TSQR has been used to accelerate the SVD step in RPCA for tall-skinny matrices AndersonTSQR , we extend this idea to the rSVD algorithm on general matrices, and we compare this communication-avoiding approach to Split-SPCP for analyzing fMRI brain scan data.

For the experiments in this section, we implement the TSQR algorithm on a single GPU card as well as multiple CPUs. The GPU card and CPUs are the same models described in Section 5.1. We use 32 CPU cores: two sockets, eight cores per socket, and two threads per core. To run QR decompositions in parallel on the GPU, we use the batched QR algorithm of MAGMA version 1.7.0 Magma , compiled as a MEX file using GCC for Red Hat version 4.8.5-11. We wrote our own implementation of TSQR on multiple CPUs using MATLAB’s Parallel Computing Toolbox. The initial number of blocks in our TSQR implementation is 32, corresponding to the 32 CPU cores we have available. Specific details on the TSQR algorithm can be found in TSQR .

fMRI brain scan data sets are inherently four-dimensional, so unwrapping the dataset into a two-dimensional array creates a tall-skinny matrix with hundreds of times more rows than columns. The enormity of these data sets makes low-rank recovery models intractable without parallelization. Figure 11 shows the total time spent performing SPCP on one of these data sets, as well as the total time spent performing rSVDs, performing the QR step of the rSVD, and moving data (labeled as the “overhead” cost). The results of this test are included in Figures 12 and 13.

The datasets used for testing were taken from a study of the human brain’s response to uncomfortable stimuli, specifically, the feeling of submerging a hand in hot or cold water. Analyzing these scans to find common neurological responses to these stimuli is difficult due to the enormous amount of error in the data. There is uniformly distributed noise due to constant ancillary physiological activity, and there are also sparsely distributed groups representing neurological structures that should all exhibit the same behavior. The ventricles, for example, are filled with cerebrospinal fluid (CSF), which does not contribute to neurological communication, so they should not be active. All signals observed in the ventricles should be treated as sparsely structured outliers. SPCP removes the uniform noise and, most remarkably, correctly identifies signals in the brain’s ventricles as outliers. In Figure 12, the largest ventricles are the two structures in the center of the brain. The rightmost image shows that the majority of the noise contained in $S$ is from these ventricles.

The other two major components of the brain are white and gray matter. The activity we are hoping to observe takes place in the gray matter, so ideally SPCP would remove most signals from the white matter regions. However, the regions of white matter are more difficult to distinguish than the regions of CSF, and SPCP removes about equal amounts of noise from the white matter as it does from the gray. If we let $S_{gm}$ be the gray-matter component of $S$, and define $S_{wm}$ and $S_{csf}$ similarly, Figure 13 shows the average BOLD signal in $S_{gm},S_{wm}$, and $S_{csf}$ for each frame in time. These data were normalized by the average original signal in the corresponding regions.

It is clear that $S_{csf}$ contains more signal than the other two regions. For $S_{gm}$, SPCP detects noise in only the first 100 time slices. The removed signal from the white matter is more distributed over time, and the total amount of noise in $S_{gm}$ and $S_{wm}$ is comparable. These results suggest that SPCP correctly identifies outliers in the fMRI data, especially within the regions of CSF.

Although using TSQR to reduce communication costs provides noticeable speedup, the non-convex approach is faster. In Figure 14, we compare our non-convex Split-SPCP algorithm to a convex solver using TSQR for the rSVD step. Both solvers were decomposing the same $106,\!496\times 462$ brain-scan data-matrix used in section 6.4, and both were run on the GPU.

As with previous tests, a reference solution $(L_{ref},S_{ref})$ was found by solving LagQN to high accuracy, and the error is measured as the normalized distance from the reference objective value. The rank of $L_{ref}$ was 423, and Split-SPCP was initialized with the rank bound $k=462$. Since $X$ had 462 columns, this was the largest possible rank bound, and Split-SPCP still greatly outperformed the convex solvers. Both solvers also found solutions of similar quality. At an error of $5\times 10^{-4}$, the low-rank component found by LagQN had rank 427, and the low-rank component found by Split-SPCP had rank 425. Similarly, the sparse component of LagQN had sparsity $58.5\%$, and the sparse component of Split-SPCP had sparsity $58.3\%$, while $S_{ref}$ was $58.5\%$ sparse. Split-SPCP recovered nearly the same solution as the parallelized convex algorithm at a much smaller cost.

6.2 Frank-Wolfe with Marginalization

Our approach to accelerating low-rank recovery models is closely connected to the one presented in FWforRPCA to accelerate the Frank-Wolfe method for RPCA. The Frank-Wolfe method, first proposed in FW , has recently gained attention for its ability to yield scalable algorithms for optimization problems involving structure-encouraging regularizers (JaggiFrankWolfe, ; FWforRPCA, ). Frank-Wolfe methods optimize a Lipschitz-differentiable function over a convex set by solving a sequence of linear subproblems. For example, if $f$ is Lipschitz-differentiable and $\mathcal{C}$ is a convex set, Frank-Wolfe solves the problem

[TABLE]

by linearizing $f$ about each iteration $x_{k}$: f(y)\approx f(x_{k})+\big{<}\nabla f(x),y-x\big{>}, minimizing the linear problem to find a decent direction $y_{k}-x_{k}$, and then stepping in this direction with a step size $\eta$: $x_{k+1}=x_{k}+\eta(y_{k}-x_{k})$. Traditionally, $\eta=\frac{2}{k+2}$ is fixed for convergence guarantees, but better performance can be achieved with adjustable step sizes FWforRPCA .

For optimization problems involving the nuclear norm, the set $\mathcal{C}$ is a scaled nuclear norm ball, and the linear subproblem becomes a search for the leading singular tuple of $L_{k}$. This is much cheaper than computing the full SVD at each iteration, which is required by proximal methods.

In contrast to nuclear norm minimization problems, it is often difficult to scale the Frank-Wolfe method for problems involving $\ell_{1}$ regularization. For these problems, each Frank-Wolfe iteration updates iterate $S_{k}$ with only a single non-zero element (FWforRPCA, ). This makes Frank-Wolfe infeasible when the sparse component of the data-matrix has several thousand non-zero elements or more. The authors of FWforRPCA adapt the Frank-Wolfe method to solve SPCP, using the traditional Frank-Wolfe update for the low-rank component and using a projected gradient update for the sparse component. Their adaptations allow for significant speedup and better convergence compared to the traditional Frank-Wolfe.

The techniques used in FWforRPCA can be seen as a special case of the marginalization we present in this work. With the sparsely structured component marginalized, the program depends only upon the low-rank component, so the benefits of the Frank-Wolfe scheme are preserved and its drawbacks are negated. Algorithm 1 shows how marginalization can be used to extend the adaptations presented in Algorithms 5 and 6 of FWforRPCA . Algorithm 1 solves the following program, which is equivalent to SPCP:

[TABLE]

With the $S$ variable marginalized, it does not affect the linear subproblems arising in the Frank-Wolfe scheme, making this approach much more scalable. Each iteration of Algorithm 1 updates the low-rank component and the nuclear norm bound, $t$, while the sparse component remains only implicitly defined. More details on this approach can be found in FWforRPCA .

While Algorithm 1 is similar to the methodology in FWforRPCA , fully marginalizing the sparse component eliminates the linear subproblem required in FWforRPCA to update the iterate $S_{k}$. We show the performance of this adapted Frank-Wolfe scheme on SPCP in Figure 15. Although its performance is better than the traditional Frank-Wolfe (see FWforRPCA ), it is still much slower than solving Split-SPCP with L-BFGS.

7 Conclusion

In this manuscript, we present a general framework for accelerating regularized low-rank recovery models on parallel architectures. Factorizing the low-rank component $L=UV^{T}$ and optimizing over the factors $U$ and $V$ induces a low-rank constraint on $L$ while eliminating the iterative SVDs often required for nuclear norm minimization. Although this factorization sacrifices the convexity of the program, we develop a certificate to determine whether our model has converged to a minimum or to a spurious stationary point, showing empirically and theoretically that spurious stationary points do not pose a serious problem for our non-convex model. We also show that marginalizing the regularized $S$ variable creates an objective that is Lipschitz-differentiable when the loss function is strongly convex. The smooth objective allows the use of first-order solvers, which exhibit faster convergence and better parallelization.

We demonstrate the effectiveness of our framework on SPCP, showing that Split-SPCP runs an order of magnitude faster on the GPU than existing solvers. There are other regularized low-rank recovery processes that decrease the cost of the SVD step for improved efficiency, including communication-avoiding and Frank-Wolfe-type algorithms. To compare to Split-SPCP, we propose an improved algorithm in both classes, and show that our SVD-free approach offers greater speedup.

With the acceleration our framework provides, regularized low-rank recovery models achieve new applications. Split-SPCP can process video streams in real-time and decompose extremely large data sets, such as fMRI brain scans, whose analysis with low-rank recovery models would otherwise be infeasible.

8 Acknowledgments

We would like to thank Tor Wager and Stephan Geuter from the Neuroscience department of the University of Colorado-Boulder for the fMRI data and advice.

Appendix A: Factorization Theorem Details

Both problems share the same set of local minimizers. This was first observed in the case of SDP in [10] and later generalized in [4].

Our variant of the theorem:

Theorem 2**.**

Consider an optimization problem of the following form

[TABLE]

where $X\in\mathbb{R}^{n\times n}$ is a positive semidefinite real matrix, and $f$ is a lower semi-continuous (lsc) function mapping to $[-\infty,\infty]$ and has a non-empty domain over the positive semi-definite matrices. Using the change of variable $X=PP^{T}$, take $P\in\mathbb{R}^{n\times k}$, and consider the problem

[TABLE]

Let $\bar{X}=\bar{P}\bar{P}^{T}$, where $\bar{X}$ is feasible for (30). Then $\bar{X}$ is a local minimizer of (30) if and only if $\bar{P}$ is a local minimizer of (31).

Proof.

We follow ideas from both [10] and [4]. From Lemma 2.1 in [10], if both $P$ and $K$ are $n\times k$ matrices, then $PP^{T}=KK^{T}$ if and only if $P=QK$ for some orthogonal matrix $Q\in\mathbb{R}^{k\times k}$. The objective in (31) depends only on $PP^{T}$, so it is clear that $P$ is a local minimizer of (31) if and only if $PQ$ is a local minimizer for all orthogonal $Q$.

Note that $g$ defined by $P\mapsto f(PP^{T})$ is also lsc since it is the composition of $f$ with the continuous function $P\mapsto PP^{T}$. We require the functions to be lsc so that the notion of “local minimizer” is well-defined.

Suppose $\bar{X}$ is a local minimizer of (30), i.e., for some $\epsilon>0$ there is no better feasible point, and factor it as $\bar{X}=\bar{P}\bar{P}^{T}$ for some $n\times k$ matrix $\bar{P}$. Then we claim $\bar{P}$ is a local minimizer of (31). If it is not, then there exists a sequence $P_{n}\rightarrow\bar{P}$ with $g(P_{n})<g(\bar{P})$ for all $n$. By continuity of the map $P\mapsto PP^{T}$, there is some $n$ large enough such that $X_{n}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}P_{n}P_{n}^{T}$ is within $\epsilon$ of $\bar{X}$, with $f(X_{n})<f(\bar{X})$, contradicting the local optimality of $\bar{X}$.

We prove the other direction using contrapositive. Suppose $\bar{X}$ is not a local minimizer of (30), so there is a sequence $X_{n}\rightarrow\bar{X}$ with $f(X_{n})<f(\bar{X})$. Factor each $X_{n}$ as $X_{n}=P_{n}P_{n}^{T}$, and observe that it is not true that $P_{n}$ converges. However, $(X_{n})$ converges and hence is bounded, thus $(P_{n})$ is bounded as well (for the spectral norm, $\|X_{n}\|=\|P_{n}\|^{2}$, and over the space of matrices, all norms are equivalent). Since $P_{n}$ are in a finite dimensional space, so the Bolzano-Weierstrass theorem guarantees that there is a sub-sequence of $(P_{n})$ that converges. Let the limit of the sub-sequence be $\bar{P}=\lim_{k\rightarrow\infty}P_{n_{k}}$, and note $\bar{P}\bar{P}^{T}=\bar{X}$ since $X_{n_{k}}\rightarrow\bar{X}$. Then $g(P_{n_{k}})=f(X_{n_{k}})<f(\bar{X})=g(\bar{P})$, so $\bar{P}$ is not a local solution. It also follows from the first paragraph of the proof that there can not be another local solution $\tilde{P}$ that also satisfies $\tilde{P}\tilde{P}^{T}=\bar{X}$. ∎

Remark 1**.**

We recover the settings of both [4, Thm. 4.1] and [10], since allowing $f$ to be extended valued and lsc encompasses constraints.

Remark 2**.**

[4, Cor. 4.2]** can be clarified to state that the two problems listed there are “equivalent” in the sense that they share the same local minimizers (as well as global), using similar continuity arguments on the mapping $\mathcal{R}$ and its adjoint. Furthermore, since the constraints are compact, solutions exist in both formulations.