Unique sparse decomposition of low rank matrices
Year: 2022

Abstract:
The problem of finding the unique low dimensional decomposition of a givenmatrix has been a fundamental and recurrent problem inmany areas. In this paper, we study the problem of seeking a unique decomposition of a low rankmatrixY ∈ R that admits a sparse representation. Specifically, we consider Y = AX ∈ R where the matrixA ∈ R has full column rank, with r < min{n, p}, and the matrixX ∈ R is element-wise sparse. We prove that this sparse decomposition ofY can be uniquely identified by recovering ground-truthA column by column, up to some intrinsic signed permutation. Our approach relies on solving a nonconvex optimization problem constrained over the unit sphere. Our geometric analysis for the nonconvex optimization landscape shows that any strict local solution is close to the ground truth solution, and can be recovered by a simple datadriven initialization followed with any second order descent algorithm. At last, we corroborate these theoretical results with numerical experiments.

heading: 1 Introduction
text: The problem of matrix decomposition has been a popular and fundamental topic under extensive investigations across several disciplines, including signal processing, machine learning, natural language processing [10, 11, 31, 46, 32, 8]. From the decomposition, one can construct efficient representation of the original data matrix. However, for any matrix Y ∈ Rp×n that can be factorized as a product of two matrices A ∈ Rp×r and X ∈ Rr×n, there exist infinitely many decompositions, simply because one can use any r× r invertible matrix Q to construct A′ = AQ and X ′ = Q−1X such that Y = AX = A′X ′, while A′ 6= A andX ′ 6= X . Thus, in various applications, additional structures and priors are being exploited to find a preferred representation [22, 15]. For example, principal component analysis (PCA) aims to find orthogonal representations which retain as much variations in Y as possible [17, 23], whereas independent component analysis (ICA) targets the representations of statistically independent non-Gaussian signals [26].
In this paper, we are interested in finding a unique sparse low-dimensional representation of Y . To this end, we study the decomposition of a low rank matrix Y ∈ Rp×n that satisfies
Y = AX, (1.1)
where A ∈ Rp×r is an unknown deterministic matrix, with r < min{n, p}, and X ∈ Rr×n is an unknown sparse matrix.
Formulation (1.1) is an important model problem in many applications. As columns of

heading: Y are viewed as linear combinations of columns of A with X being the sparse coefficient, (1.1) can be used to form overlapping clusters of the n columns of Y via the support of
text: X with columns of A being viewed as r cluster centers [12, 7]. When we form a p × r low-dimensional representation of Y via sparse combinations, this greatly enhance the interpretability of the resulting representations [28, 21, 4], in the same spirit as the sparse PCA, but (1.1) generalizes to the factorization of non-orthogonal matrices.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
To motivate our approach, we first consider the simple case that A has orthonormal columns, namely, ATA = Ir1. Then it is easy to see that the sparse coefficient matrix X is recovered by multiplying Y on the left byAT ,
ATY = ATAX = X. (1.2)
The problem of finding such orthonormal matrix A boils down to successively finding a unit-norm direction q that renders qTY as sparse as possible [34, 39, 35],
min q
∥ ∥qTY ∥ ∥
sparsity s. t. ‖q‖2 = 1. (1.3)
However, the natural choice of sparsity penalty, either ℓ0 or ℓ1, leads to trivial and meaningless solutions, as there always exists q in the null space ofAT such that qTY = 0.
To avoid the null space ofAT , we instead choose to find the unit direction q that maximizes the ℓ4 norm of qTY as
max q
∥ ∥qTY ∥ ∥
4 s. t. ‖q‖2 = 1. (1.4)
The above formulation is based on the key observation that the objective value ismaximized when q coincideswith one column ofA (see, Section 2, for details) while the objective value is zero when q lies in the null space of AT . The ℓ4 norm objective function and its variants have been adopted as a sparsity regularizer in a line of recent works [30, 44, 43, 35, 42]. However, even with this new objective function, the null space ofAT persists as a challenge for solving the optimization problem: they form a flat region of saddle points.
This paper characterizes the nonconvex optimization landscape of (1.4) and proposes a guaranteed procedure that avoids the flat region and provably recovers the global solution to (1.4), which corresponds to one column of A. More specifically, we demonstrate that, despite the non-convexity, (1.4) still possesses benign geometric property in the sense that any strict local solution with large objective value is globally optimal and recovers one column ofA, up to its sign. See, Theorem 3.1 in Section 3.1 for the population level result and Theorem 3.4 for the finite sample result.
We further extend these results to the general case when A only has full column rank in Theorem 3.6 of Section 3.2. To recover a general A with full column rank, our procedure first resorts to a preconditioning procedure of Y proposed in Section 2.3 and then solves a optimization problem similar to (1.4). From our analysis of the optimization landscape, the intriguing problem boils down to developing algorithms to recover the nontrivial local solutions by avoiding regions with small objective values. We thus propose a simple initialization scheme in Section 4.1 and prove in Theorem 4.3 that such initialization, proceeded with any second order descent algorithm [20, 27], suffices to find the global solution, up to some statistical error. Our theoretical analysis provides the explicit convergence rate of the statistical error and characterizes its dependence on various dimensions, such as p, r and n, as well as the sparsity of X .
Numerical simulation results are provided in Section 5. Due to the space limitation, we defer all the proof along with our conclusions and discussion of several future directions of our work to Appendix.
Notations Throughout this paper, we use bold lowercase letters, like a, to represent vectors and bold uppercase letters, like A, to represent matrices. For matrix X , Xij denotes the entry at the i-th row and j-th column ofX , withXi· andX·j denoting the i-th row and j-th column ofX , respectively. Oftentimes, we writeX·j = Xj for simplicity. We use grad and Hess to represent the Riemannian gradient and Hessian. For any vector v ∈ Rd, we use ‖v‖q to denote its ℓq norm, for 1 ≤ q ≤ ∞. The notation v◦q stands for {v q i }i. For matrices, we use ‖·‖F and ‖·‖op to denote the Frobenius norm and the operator norm, respectively. For any positive integer d, wewrite [d] = {1, 2, . . . , d}. The unit sphere in d-dimensional real space Rd is written as Sd−1. For two sequences an and bn, we write an À bn if there exists some constant C > 0 such that an ≤ Cbn for all n. Both uppercase C and lowercase c are reserved to represent numerical constants, whose values may vary line by line.
1 Ir is the identity matrix of size r × r.

heading: 1.1 Related work
text: Finding the unique factorization of a matrix is an ill-posed problem in general due to infinitely many solutions. There exist several strands of studies from different contexts on finding the unique decomposition of Y by imposing additional structures onA andX . We start by reviewing the literature which targets the sparse decomposition of Y .
Dictionary learning The problems of dictionary learning (DL) [2, 38, 19, 39] and sparse blind deconvolution or convolutional dictionary learning [13, 29] study the unique decomposition of Y = AX whereX is sparse andA has full row rank. In this case, the row space of Y lies in the row space ofX , suggesting to recover the sparse rows ofX via solving the following problem,
min q
∥ ∥qTY ∥ ∥
1 s. t. q 6= 0. (1.5)
Under certain scaling and incoherence conditions onA, the objective achieves theminimum value when q is equal to one column ofA, at the same time qTY recovers one sparse row of X . This idea has been studied andmodified in a strand of papers whenA has full row rank [38, 39, 45, 30, 44, 35, 42, 37, 47]. In our context, the major difference rises in the matrix A, which has full column rank rather than row rank, therefore minimizing ∥ ∥qTY ∥ ∥
1 as before
only leads to some vector in the null space of AT , yielding the trivial zero objective value.
We would love to note that [35] uses the same objective function in (1.4) to study the problem of overcomplete dictionary learning (where A has full row rank), however the optimization landscape when A has full column rank is significantly different from that in the overcomplete setting. The more complicated optimization landscape in our setting brings additional difficulty of the analysis and requires a proper initialization in our proposed algorithm. We refer to Appendix B for detailed technical comparison with [35].
Sparse PCA Sparse principal component analysis (SPCA) is a popular method that recovers a unique decomposition of a low-rank matrix Y by utilizing the sparsity of its singular vectors. However, as being said, under Y = AX , SPCA is only applicable when X coincides with the right singular vectors of Y . Indeed, one formulation of SPCA is to solve
max U∈Rn×r
tr ` UTY TY U ˘ − λ‖U‖1, s. t. U TU = Ir, (1.6)
which is promising only if X corresponds to the right singular vectors of Y . It is worth mentioning that among the various approaches of SPCA, the following one might be used to recover one sparse row of X ,
min u,v
∥ ∥Y − uvT ∥ ∥
2 2 + λ‖v‖1 s. t. ‖u‖2= 1. (1.7)
This procedure was originally proposed by [49] and [36] together with an efficient algorithm by alternating the minimization between u and v. However, there is no guarantee that the resulting solution recovers the ground truth.
Factor analysis Factor analysis is a popular statistical tool for constructing low-rank representations of Y by postulating Y = AX + E where A ∈ Rp×r is the so-called loading matrix with r = rank(A) < min{n, p},X ∈ Rr×n contains n realizations of a r-dimensional factor and E is some additive noise. Here only Y is observable. Factor analysis is mainly used to recover the low-dimension column space ofA or the row space ofX , rather than to identify and recover the unique decomposition. Recently, [7] studied the unique decomposition of Y when the columns of X are i.i.d. realizations of a r-dimensional latent random factor. The unique decomposition is further used for (overlapping) clustering the rows of Y via the assignment matrix A. To uniquely identify A, [7] assumes that A contains at least one r × r identity matrix, coupled with other scaling conditions onA (we refer to [7] for detailed discussions of other existing conditions in the literature of factor models that ensure the unique decomposition of Y but require strong prior information on either A or X). By contrast, we rely on the sparsity of X instead of A which is more general than requiring the existence of a r × r identity matrix.
NMF and topic models Such existence condition of identity matrix in either A or X has a variant in non-negative matrix factorization (NMF) [14] and topic models [3, 8, 9], also
see the references therein, where Y , A and X have non-negative entries. Since all Y , A andX from model (1.1) are allowed to have arbitrary signs in our context, the approaches designed for NMF and topic models are inapplicable.

heading: 2 Formulation and Assumptions
text: The decomposition of Y = AX is not unique without further assumptions. To ensure the uniqueness of such decomposition, we rely on two assumptions on the matrices A and X , stated in the following Section 2.1.
Our goal is to uniquely recoverA fromY , up to a some signed permutation. More precisely, we aim to recover columns of AP for some signed permutation matrix P ∈ Rr×r. To facilitate the understanding and motivate our approach, in Section 2.2 we first state our procedure for the unique recovery of A when A has orthonormal columns. Its theoretical analysis is presented in Section 3. Later in Section 2.3, we discuss how to extend our results to the case when A is a more general full column rank matrix under Assumption 2.2.
For now, we only focus on the recovery of one column ofA as the remaining columns can be recovered via the same procedure after projecting Y onto the complement space spanned by the recovered columns of A (see Section D for detailed discussion).

heading: 2.1 Assumptions
text: We first resort to the matrix X ∈ Rr×n being element-wise sparse. The sparsity of X is modeled via the Bernoulli-Gaussian distribution, stated in the following assumption.
Assumption 2.1 Assume Xij = BijZij for i ∈ [r] and j ∈ [n], where
Bij i.i.d. ∼ Ber(θ), Zij i.i.d. ∼ N (0, σ2). (2.1)
The Bernoulli-Gaussian distribution is popular for modeling sparse random matrices [38, 2, 1, 39]. The overall sparsity level of X is controlled by θ, the parameter of the Bernoulli distribution. We remark that the Gaussianity is assumed only to simplify the proof and to obtain more transparent deviation inequalities between quantities related withX and their population counterparts. Both our approach and analysis can be generalized to cases where Zij are centered i.i.d. sub-Gaussian random variables.
We also need another condition on the matrix A. To see this, note that even when A were known, recovering X from Y = AX requires A to have full column rank. We state this in the following assumption.
Assumption 2.2 Assume the matrixA ∈ Rp×r has rank(A) = r with ‖A‖op= 1.
The unit operator normofA is assumedwithout loss of generality as one can always re-scale σ2, the variance of X , by ‖A‖op.

heading: 2.2 Recovery of the orthonormal columns of A
text: In this section, we consider the recovery of one column of A when A is a semi-orthogonal matrix satisfying the following assumption.
Assumption 2.3 Assume ATA = Ir.
Our approach recovers columns of A one at a time by adopting the ℓ4 maximization to penalize the sparsity of rows of matrix X . Its rationale is based on the following lemma, assuming the orthogonality among columns of A.
Lemma 2.4 Under Assumption 2.3, solving the following problem
max q
∥ ∥ATq ∥ ∥
4 4 s. t. ‖q‖2 = 1 (2.2)
recovers one column of A, up to its sign.
Intuitively, under Assumption 2.3, we have ‖ATq‖2≤ 1 for any unit vector q. Therefore, criterion (2.2) seeks a vector ATq within the unit ball to maximize its ℓ4 norm. When q corresponds to one column of A, that is, q = ai for any i ∈ [r], we have the largest objective ‖ATai‖44= 1. This ℓ4 norm maximization approach has been used in several related literature, for instance, sparse blind deconvolution [44, 30], complete and over-complete dictionary learning [43, 42, 35], independent component analysis [25, 24] and tensor decomposition [18].
The appealing property of maximizing the ℓ4 norm is its benign geometry landscape under the unit sphere constraint. Indeed, despite of the non-convexity of (2.2), our result in Theorem 3.1 implies that any strict location solution to (2.2) is globally optimal. This enables us to use any second order gradient ascent method to solve (2.2). Motivated by Lemma 2.4, since we only have access to Y ∈ Rp×n, we propose to solve the following problem to recover one column of A,
min q
F (q) . = −
1
12θσ4n
∥ ∥Y Tq ∥ ∥
4 4 s. t. ‖q‖2 = 1. (2.3)
The scalar (12θσ4n)−1 is a normalization constant. The following lemma justifies the usage of (2.3) and also highlights the role of the sparsity of X .
Lemma 2.5 Under model (1.1) and Assumption 2.1, we have
E rF (q)s = −1 4
” (1− θ) ∥ ∥ATq ∥ ∥ 4
4 + θ
∥ ∥ATq ∥ ∥
4
2
ı (2.4)
where the expectation is taken over the randomness of X .
Remark 2.6 (Role of the sparsity parameter θ) Lemma 2.5 implies that, for large n, solving (2.3) approximately finds the solution to
min q f pqq .= −1 4
” (1− θ) ∥ ∥ATq ∥ ∥ 4
4 + θ
∥ ∥ATq ∥ ∥
4
2
ı s. t. ‖q‖2 = 1 (2.5)
The objective function is a convex combination of ‖ATq‖44 and ‖A Tq‖22 with coefficients depending on the magnitude of θ. In view of Lemma 2.4, it is easy to see that solving (2.5) recovers one column of A, up to the sign, as long as θ < 1. However, the magnitude of θ controls the benignness of the geometry landscape of (2.5). When θ is small, orX is sufficiently sparse, we essentially solve (2.2) which has the most benign landscape. On the other hand, when θ → 1, the landscape of (2.5) is mostly determined by the eigenvalue problem2 which maximizes ‖ATq‖2 subject to ‖q‖2= 1. We will demonstrate that when X is sufficiently sparse, second order descent algorithm with a simple initialization finds the globally optimal solution to (2.3) in Section 3.

heading: 2.3 Recovery of the non-orthogonal columns of A
text: In this section, we discuss how to extend our procedure to recoverA from Y = AX when A is a general full column rank matrix satisfying Assumption 2.2. The main idea is to first resort to a preconditioning procedure of Y such that the preconditioned Y has the decomposition ĀX̄ , up to some small perturbation, where Ā satisfies Assumption 2.3 and X̄ satisfies Assumption 2.1 with σ2 = 1. Then we apply our procedure in Section 2.2 to recover Ā. The recovered Ā is further used to recover the original A.
To precondition Y , we propose to left multiply Y by the following matrix
D . = ”` Y Y T ˘+ı1/2 ∈ Rp×p (2.6)
where M+ denotes the Moore-Penrose inverse of any matrix M . The resulting preconditioned Y satisfies
Ȳ . = DY = ĀX̄ +E (2.7)
with Ā satisfying Assumption 2.3, X̄ = X/ ? θnσ2 andE being a perturbation matrix with small entries. We refer to Proposition 3.5 below for its precise statements.
Analogous to (2.3), we propose to recover one column of Ā by solving the following problem
min ‖q‖
2 =1
Fg(q) . = −
θn
12
∥ ∥Ȳ Tq ∥ ∥
4 4 (2.8)
Theoretical guarantees of this procedure are provided in Section 3.2. After recovering one column of Ā, the remaining columns of Ā can be successively recovered via the procedure in Section D. In the end,A can be recovered by first inverting the preconditioning matrixD as D−1Ā and then re-scaling its largest singular value to 1.
2When A is orthonormal, this eigenvalue problem processes the worst landscape as there are infinitely many solutions obtaining the same eigenvalue.

heading: 3 Theoretical Guarantees
text: We provide theoretical guarantees for our procedure (2.3) in Section 3.1 when A has orthonormal columns. The theoretical guarantees of (2.8) for recovering a general full column rank A are stated in Section 3.2.

heading: 3.1 Theoretical guarantees for semi-orthonormal A
text: In this section, we provide guarantees for our procedure by characterizing the solution to (2.3) when A satisfies Assumption 2.3.
As the objective function F (q) in (2.3) concentrates around f(q) in (2.5), it is informative to first analyze the solution to (2.5). Although (2.5) is a nonconvex problem and has multiple local solutions, Theorem 3.1 below guarantees that any strict local solution to (2.5) is globally optimal, in the sense that, it recovers one column ofA, up to its sign. We introduce the null region R0 of our objective in (2.5),
R0 = { q ∈ Sp−1 : ‖ATq‖∞= 0 } . (3.1)
Theorem 3.1 (Population case) Under Assumption 2.3, assume θ ≤ 1/6. Any local solution q̄ to (2.5), that is not in R0, satisfies
ā = APe1 (3.2)
for some signed permutation matrix P ∈ Rr×r.
The detailed proof of Theorem 3.1 is deferred to Appendix F.3. We only offer an outline of our analysis below.
The proof of Theorem 3.1 relies on analysis of the optimization landscape of (2.5) on disjoint partitions of Sp−1 = {q ∈ Rp : ‖q‖2= 1}3, defined as
R1 . = R1(C⋆) = { q ∈ Sp−1 : ‖ATq‖2∞≥ C⋆ } , (3.3) R2 = S p−1 \ pR0 ∪R1q .
HereC⋆ is any fixed constant between 0 and 1. The upper bound follows from the inequality that ‖ATq‖∞= maxk|aTk q|≤ ‖ak‖2‖q‖2= 1 for any q ∈ S
p−1. The region R0 can be easily avoided by choosing the initialization such that the objective function f(q) is not equal to zero. For R1 and R2, we are able to show the following results. Let Hess f(q) be the Riemannian Hessian matrix of (2.5) at any point q ∈ Sp−1.
(1) Optimization landscape for R1:
Lemma 3.2 Assume θ < 1. Any local solution q̄ ∈ R1(C⋆) to (2.5) with C⋆ > 1 2
b θ
1−θ
recovers one column of A, that is, for some signed permutation matrix P
q̄ = AP·1.
Lemma 3.2 shows that any critical point q ∈ R1 is either a strict saddle point that there exists a direction along which the Hessian is negative, or the desired local solution q̄ that satisfies the second order optimality condition and is equal to one column of A, up to its sign.
(2) Optimization landscape for R2:
Lemma 3.3 Assume θ < 1/3. For any point q ∈ R2(C⋆) with C⋆ ≤ 1−3θ
2 , there exists v such that
vT Hess f pqqv < 0. (3.4)
Lemma 3.3 implies that any critical point inR2 is a saddle point that can be escaped by negative curvature. Hence there is no local solution to (2.5) in the region R2.
3Visualization of the partitions in S2 is available in section A.
Theorem 3.1 thus follows from Lemma 3.2 and Lemma 3.3, provided that c
θ
1− θ < 1− 3θ. (3.5)
Condition (3.5) puts restrictions on the upper bound of θ. It is easy to see that (3.5) holds for any θ ≤ 1/6. As discussed in Remark 2.6, a smaller θ leads to amore benign optimization landscape. In light of Theorem 3.1, we now provide guarantees for the solution to the finite sample problem (2.3) in the following theorem. Define the sample analogue of the null region R0 in (3.1) as
R′0(c⋆) . = { q ∈ Sp−1 : ‖ATq‖2∞ ≤ c⋆ }
(3.6)
for any given value c⋆ ∈ [0, 1). Theorem 3.4 (Finite sample case) Under Assumptions 2.1 and 2.3, assume θ ∈ (0, 1/9] and
n ≥ Cmax
{
r2 c⋆ , log2 n
}
r log n
θc⋆ (3.7)
for some sufficiently large constant C > 0 and any c⋆ ∈ (0, 1/4]. Then with probability at least 1− cn−c ′
, any local solution q̄ to (2.3) that is not in R′0(c⋆) satisfies
‖q̄ −AP·1‖ 2 2 À
c r2 log n
θn +
ˆ θr2 + log2 n
θ
˙ r log n
n (3.8)
for some signed permutation matrix P . Here we defer our discussion of technical details and full proof in section C.

heading: 3.2 Theoretical guarantees for general full column rank A
text: In this section, we provide theoretical guarantees for our procedure of recovering a general full column rank matrix A under Assumption 2.2.
Recall from Section 2.3 that our approach first preconditions Y by usingD from (2.6). The following proposition provides guarantees for the preconditioned Y , denoted as Ȳ = DY . The proof is deferred to Appendix F.5. Write the SVD of A = UADAV TA with UA ∈ R p×r and VA ∈ Rr×r being, respectively, the left and right singular vectors. Proposition 3.5 Under Assumptions 2.1 and 2.2, assume n ≥ Cr/θ2 for some sufficiently large constant C > 0. With probability greater than 1− 2e−c ′r, one has
Ȳ = ĀX̄ +E (3.9)
where Ā = UAV T A , X̄ = X/
? θnσ2 and E = Ā∆X̄ with
‖∆‖op ≤ c ′′ 1
θ
c r
n . (3.10)
Here c′ and c′′ are positive constants. Proposition 3.5 implies that, when n ≥ Cr/θ2, the preconditioned Y satisfies
Ȳ = Ā(Ir +∆)X̄ ≈ ĀX̄ (3.11)
with ĀT Ā = Ir. This naturally leads us to apply our procedure in Section 2.2 to recover columns of Ā via (2.8). We formally show in Theorem 3.6 below that any local solution to (2.8) approximately recover one column of Ā up to a signed permutation matrix. Similar to (3.6), define
R ′′ 0 (c⋆) . =
{
q ∈ Sp−1 : ∥ ∥ĀTq ∥ ∥
2 ∞ ≤ c⋆
}
(3.12)
for some given value c⋆ ∈ [0, 1). Theorem 3.6 Under Assumption 2.1 and 2.2, assume θ ∈ p0, 1/9s and
n ≥ C r
c⋆θ max
{
log3 n, log n c⋆θ ? θ , log2 n c⋆θ , r c⋆ ? θ , r2 log n c⋆ } . (3.13)
Then with probability at least 1 − cn−c ′ − 4e−c ′′r, any solution q̄ to (2.8) that is not in Region R ′′
0 pc⋆q satisfies ∥
∥q̄ − ĀP·1 ∥ ∥ 2 2 À
c r log n
θ2n +
c r2 log n
θn +
ˆ θr2 + log2 n
θ
˙ r log n
n (3.14)
for some signed permutation matrix P .
The proof of Theorem 3.6 can be found in Appendix F.6. Due to the preconditioning step, the requirement of the sample size in (3.13) is slightly stronger than (3.7), whereas the estimation error of q̄ only has an additional a r log n/(θ2n) term comparing to (3.8).
Theorem 3.6 requires to avoid the null region R′′0 (c⋆) in (3.12). We provide a simple initialization in the next section that provably avoids R ′′
0 . Furthermore, every iterate of any
descent algorithm based on such initialization is provably not in R ′′
0 either.

heading: 4 Complete Algorithm and Provable Recovery
text: In this section, we present a complete pipeline for recovering A from Y . So far we have established that every local solution to (2.8), that is not in R′′0 (c⋆), approximately recovers one column of Ā = UAV TA . To our end, we will discuss: (1) a data-driven initialization in Section 4.1 which, together with Theorem 3.6, provably recovers one column of Ā; (2) a deflation procedure [38, 39, 35] in Section D that sequentially recovers all remaining columns of Ā. Due to the limitation of space we defer our discussion of deflation procedure in appendix.

heading: 4.1 Initialization
text: Our goal is to provide a simple initialization such that solving (2.8) via any second order descent algorithm provably recovers one column of Ā. According to Theorem 3.6, such an initialization needs to guarantee the following conditions.
• Condition I: The initial point q(0) does not fall into region R ′′
0 (c⋆) for some c⋆ satisfying (3.13) in Theorem 3.6.
• Condition II: The updated iterates q(k), for all k ≥ 1, stay away from R′′0 (c⋆) as well.
We propose the following initialization
q(0) = Ȳ 1n
‖Ȳ 1n‖2 ∈ Sp−1. (4.1)
The following two lemmas guarantee that both Condition I and Condition II are met for this choice. Their proofs can be found in Appendices F.7 and F.8.
Lemma 4.1 Under Assumption 2.1 and 2.2, assume θ ∈ p0, 1/9s and
n ≥ C r2
θ max
{
log3 n, r log n
θ ? θ
, r log2 n
θ , r2? θ , r3 log n
}
. (4.2)
holds, then, with probability at least 1−2e−cr, the initialization q(0) in (4.1) is not in regionR ′′
0 (c⋆) with c⋆ = 1/(2r).
Lemma 4.2 Let q(k), for k ≥ 1, be any updated iterate from solving (2.3) by using any monotonic decreasing algorithm with the initial point q(0) chosen as (4.1). If
n ≥ C r2
θ max
{
log3 n, r log n
θ ? θ , r2? θ , θ2r2 log n
}
(4.3)
holds, then, with probability at least 1− cn−c ′ − 2e−c ′′r, one has
q(k) /∈ R′′0 (c⋆), for all k ≥ 1,
with c⋆ = 1/(2r).
Combining Lemmas 4.1 and 4.2 together with Theorem 3.6 readily yields the following theorem. Theorem 4.3 Under Assumptions 2.1 and 2.2, assume θ ∈ (0, 1/9] and (4.2) holds. Let q̄ be any local solution to (2.8) from any monotonic decreasing second order algorithms with the initial point chosen as (4.1). With probability at least 1− cn−c ′ − 4e−c ′′r, one has
‖q̄ − ĀP·1‖ 2 2 À
c r log n
θ2n +
c r2 log n
θn +
r log3 n
θn for some signed permutation matrix P .
Theorem 4.3 provides the guarantees for using any monotonic decreasing second order algorithms [33, 5] to solve (2.8) with the initialization chosen in (4.1).

heading: 5 Experiments
text: In this section we verify the empirical performance of our proposed procedure for recovering A under model (1.1) in different scenarios. Due to the space limit, we defer more experiments to the Appendix of this paper.
Experiment setup
To generate the data Y = AX , we generate the columns of A by using the normalized left singular vectors of R ∈ Rp×r where Rij i.i.d. ∼ N (0, 1). The sparse coefficient matrix X ∈ Rr×n are generated as Xij i.i.d. ∼ BG(θ). To evaluate the success of recovering one column vector of A for any estimate q ∈ Sp−1, we use the following criterion,
Err pqq = min 1≤i≤r p1− |〈q,ai〉 |q (5.1)
If Err pqq ≤ ρe, we say the vector q recovers the ground-truth column vector of A. We choose ρe = 1 × 10−2 in our simulation settings. To evaluate the recovery of the whole matrixA, we use the following normalized Frobenius norm between any estimateAest and the true A:
min P 1? r ‖Aest −AP ‖F s.t. P is a signed permutation matrix. (5.2)
Wefirst evaluate the probability of successfully recovering one columnofA in two scenarios. In the first case, we vary simultaneously θ and r while in the second case we change n and r. We then evaluate the performance of our procedure, Algorithm 1 in Section D, for recovering the full matrix A.
Recovery probability with varying θ and r
We fix p = 100 and n = 5×103 while vary θ ∈ {0.01, 0.04, . . . , 0.58} and r ∈ {10, 30, . . . , 70}. For each pair of pθ, rq, we repeatedly generate 200 data sets and apply our procedure in (2.8). The averaged recovery probability of our procedure over the 200 replicates is shown in Figure 1a. The recovery probability gets larger as r decreases, in line with Theorem 3.6. We also note that the recovery increases for smaller θ. This is because smaller θ renders a nicer geometric landscape of the proposed non-convex problem, as detailed in Remark 2.6. On the other hand, the recovery probability decreases when θ is approaching to 0. As suggested by Theorem 3.4, the statistical error of estimating A gets inflated as θ gets too small.
Recovery probability with varying n and r
Here we fix p = 100 and the sparsity parameter θ = 0.1. We vary r ∈ {10, 30, . . . , 70} and n ∈ {2000, 3000, . . . , 12000}. Figure 1b shows the averaged recovery probability of our procedure over 200 replicates in each setting. Our procedure performs increasingly better as n increases, as expected from Theorem 3.4.

heading: 6 Conclusion and Future Work
text: In this paper, we have studied the unique decomposition of a low rank matrix Y that admits a sparse low-dimensional representation. Under model Y = AX where X has i.i.d. Bernoulli-Gaussian entries and A has full column rank, we propose a nonconvex procedure that provably recovers A, a quantity that can be further used to recover X . We provide a complete analysis for recovering one column of A, up to the sign, by showing that any second order descent algorithm provably attains the global solution with a simple and data-driven initialization, despite the nonconvex nature of the proposed procedure.
There are several directions that are certainly worth further pursuing. For instance, a complete analysis of the deflation procedure for recovering the full matrixA is certainly of great interest. It is also worth studying this decomposition problem in presence of some additive errors, that is, Y = AX+E. Our current procedure only toleratesE that has small entries. How to modify our procedure to accommodate a moderate / large E is an interesting and challenging problem that we leave to future research.
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0
0.2
0.4
0.6
0.8
1
1.2
R e c o v e ry
p ro
b a b ili
ty
r=10 r=30 r=50 r=70
(a) Recovery probability versus θ: the averaged probability of successful recovery for different θ and rwith p = 100 and n = 1.2×104.
2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000
n
0
0.2
0.4
0.6
0.8
1
1.2
R e c o v e ry
p ro
b a b ili
ty
r=10 r=30 r=50 r=70
(b) Recovery probability versus n: the averaged probability of successful recovery for different n and r with p = 100 and θ = 0.1.

Review summaries:
This paper studies the unique decomposition of a low-rank matrix into sparse and full column-rank components. Explicit conditions under which such decomposition is possible (up to a signed permutation) are derived. Moreover, an initialization scheme is proposed that, together a generic second order method, guarantees the exact and unique recovery of the sparse and full column-rank components. The authors first study the setting where the full column-rank component has orthogonal columns. Then they extend their results to the general setting where this component is no longer column-orthogonal. To the best of my knowledge, the technical contributions of the paper are sound and correct.

The authors propose a new method and analysis to decompose a given matrix $Y$ as $Y=AX$ with $X$ sparse and $A$ full column rank. This is similar but different to sparse PCA and dictionary learning, and has interesting applications.  The proposed algorithm essentially recovers $A$ one column at a time (after deflation) and maximizes the $\ell_4$ norm of the overlap of a unit-norm vector with $Y$. As the authors also state, using the $\ell_4$ norm is not entirely new but it gives a benign optimization landscape (in particular, strict local minimizers that are close to the desired solutions).  The best results (theoretically and experimentally) are obtained when $A$ is an orthonormal matrix and $X$ is a random matrix from a Bernouilli-Gaussian distribution. The authors propose a procedure (called preconditioning) that generalizes their results (again theoretically and experimentally) to any $A$ but it is difficult to see how good the results are. In any case, the provided experiment clearly show that we do no longer have exact recovery in this case.  While the theoretical results seem -- I did not check the proofs in the 40(!) page appendix -- nontrivial and interesting, they only apply when recovering one vector and for a random $X$ from the Bernouilli-Gaussian distribution. No results are given for general $X$. In addition, the proposed deflation strategy is a typical technique to extend to multiple vector but typically performs quite poorly when $A$ is ill-conditioned (near linear depednent columns). The authors should therefore restate their contributions in the abstract and introduction to better reflect their main contributions.

The paper studies the problem of computing a low rank decomposition of a matrix. Specifically, let $A$ be a $p \times r$ matrix and $B$ be a $X \times n$ matrix. Given the matrix $Y = AX$, the authors study the problem of computing the matrices $A$, $X$. Without further assumptions, the problem is ill-defined as for any decomposition $A$, $X$ and any invertible matrix $Q$, $AQ$ and $Q^{-1}X$ is another decomposition of the matrix $Y$.   The authors make an assumption here that the matrix $X$ is sparse in particular the matrix $X$ has iid entries where each entry is independently nonzero with probability $\theta$ and conditioned on an entry being nonzero it has gaussian distribution. Assume that $A$ has orthonormal columns. The authors argue that solving the maximization problem $\max_{q: \|q\|_2 = 1}\|q^T Y\|_4^4$ recovers a column of the matrix $A$ upto a sign with some error additive error. Then they show that for the general case of $A$, they can use a preconditioner computed using $Y$ that would reduce the problem to the case  orthonormal matrix and additive noise matrix of small norm.   They argue that any 2nd order optimization algorithm can be used by arguing that all critical points are either close to recovering the optimal column or that they are saddle points and that there is a direction in which the hessian is negative.  They show some experiments by generating a random matrix $A$ with orthonormal columns and generating a matrix $X$ from the above described Gaussian-Bernoulli distribution. They show that as the matrix $X$ gets sparser, the probability of recovery of $A$ increases as predicted by their theorems and also that as $n$, the number of columns of the matrix $X$ increases, the probability of recovery of $A$ increases.

The paper studies the problem of decomposing a low-rank p x n matrix Y into Y = AX, where A is p x r of full column rank and X is a sparse r x n matrix.  The paper first considers the case where A has orthonormal columns. It shows that by accessing Y only, one can recover approximately all columns of A up to the signs. It constructs a minimization programme involving the Schatten-4 norm, whose solution is an approximation to some column of A. The minimization problem can be solved by second-order methods. Then one can project A away from this solution and repeat the process. The paper gives quantitative error guarantees, provided that n is sufficiently large. The paper generalizes the result of the orthonormal case to a general A by preconditioning.