In several problems of contemporary interest, arising for instance, in recommender system applications, for example, the Netflix Prize competition (SIGKDD and Netflix, 2007), observed data is in the form of a large sparse matrix, Ω, where Ω . Popularly dubbed as the matrix completion problem (Cand`es and Recht, 2009; Mazumder et al., 2010), the task is to predict the unobserved entries, under the assumption that the underlying matrix is of low-rank. This leads to the natural rank regularized optimization problem:

where, ) denotes the projection of onto the observed indices Ω and is zero otherwise; and denotes the usual Frobenius norm of a matrix. Problem (1), however, is computationally difficult due to the presence of the combinatorial rank constraint (Chistov and Grigor’ev, 1984). A natural convexification (Fazel, 2002; Recht et al., 2010) of rank(, the nuclear norm of X, which leads to the following surrogate of Problem (1):

Cand`es and Recht (2009); Cand`es and Plan (2010) show that under some assumptions on the underlying “population” matrix, a solution to Problem (2) approximates a solution to Problem (1) reasonably well. The estimator obtained from Problem (2) works quite well: the nuclear norm shrinks the singular values and simultaneously sets many of the singular values to zero, thereby encouraging low-rank solutions. It is thus not surprising that Problem (2) has enjoyed a significant amount of attention in the wider statistical community over the last decade. There have been impressive advances in understanding its statistical properties (Cand`es and Plan, 2010; Cand`es and Tao, 2010; Recht et al., 2010; Recht, 2011; Gross, 2011; Rohde and Tsybakov, 2011; Koltchinskii et al., 2011; Negahban and Wainwright, 2011; Chen, 2015; Lecu´e and Mendelson, 2018; Chen et al., 2019b). Motivated by the work of Cand`es and Recht (2009); Cai et al. (2010), the authors in Mazumder et al. (2010) proposed Soft-Impute, an EM-flavored (Dempster et al., 1977) algorithm for optimizing Problem (2). For some other computational work in developing scalable algorithms for Problem (2), see the papers Jaggi and Sulovsk´y (2010); Freund et al. (2015); Hastie et al. (2016), and references therein. Typical assumptions under which the nuclear norm works as a good proxy for the low-rank problem require the entries of the singular vectors of the “true” low-rank matrix to be sufficiently spread, and the missing pattern to be roughly uniform. The proportion of observed entries needs to be sufficiently larger than the number of parameters of the matrix O ((m + n)r), where, r denotes the rank of the true underlying matrix. Some extensions under general sampling distribution has been made in Klopp (2014); Alquier (2015). Negahban and Wainwright (2012) proposes improvements with a (convex) weighted nuclear norm penalty in addition to spikiness constraints for the noisy matrix completion problem.

The nuclear norm penalization framework, however, has limitations. If some conditions mentioned above fail, Problem (2) may fall short of delivering reliable low-rank estimators with good prediction performance (on the missing entries). Since the nuclear norm shrinks the singular values, in order to obtain an estimator with good explanatory power, it often results in a matrix estimator with high numerical rank — thereby leading to models that have higher rank than what might be desirable. The limitations mentioned above, however, should not come as a surprise to an expert — especially, if one draws a parallel connection to the Lasso (Tibshirani, 1996), a popular sparsity inducing shrinkage mechanism effectively used in the context of sparse linear modeling and regression. In the linear regression context, the Lasso often leads to dense models and suffers when the features are highly correlated — the limitations of the Lasso are quite well known in the statistics literature, and there have been major strides in moving beyond the convex -penalty to more aggressive forms of nonconvex penalties (Fan and Li, 2001; Zou and Li, 2008; Mazumder et al., 2011; Zhang, 2010; Zhang and Zhang, 2012; Loh and Wainwright, 2015; Bertsimas et al., 2016; Zheng et al., 2017; Feng and Zhang, 2017). The key principle in these methods is the use of nonconvex regularizers that better approximate the -penalty, leading to possibly nonconvex estimation problems. Thusly motivated, we study herein, the following family of nonconvex regularized estimators for the task of (noisy) matrix completion:

where, 1 are the singular values of ) is a concave penalty function on [0) that takes the value 0. We will denote an estimator obtained from Problem (3) by ˆThe family of penalty functions ) is indexed by the parameters (parameters together control the amount of nonconvexity and shrinkage — see for example Mazumder et al. (2011); Zhang and Zhang (2012) and also Section 2, herein, for examples of such nonconvex families.

A caveat in considering problems of the form (3) is that they lead to nonconvex optimization problems and thus obtaining a certifiably optimal global minimizer is generally difficult. Fairly recently, Bertsimas et al. (2016); Mazumder and Radchenko (2015) have shown that subset selection problems in sparse linear regression can be computed using advances in mixed integer quadratic optimization. Such global optimization methods, however, do not apply to matrix variate problems involving spectralpenalties, as in Problems (1) or (3). The main focus in our work herein is to develop a computationally scalable algorithmic framework that allows us to obtain high quality stationary points or upper boundsfor Problem (3) — we obtain a path of solutions ˆacross a grid of values of (Problem (3) by employing warm-starts, following the path-following scheme proposed in Mazumder et al. (2011). Leveraging problem structure, modern advances in computationally scalable low-rank SVDs and appropriately advancing the tricks successfully employed in Mazumder et al. (2010); Hastie et al. (2016), we empirically demonstrate the computational scalability of our method for problems of the size of the Netflix dataset, a matrix of size (approx.) 480observed entries. Perhaps most importantly, we demonstrate empirically that the resultant estimators lead to better statistical properties (i.e., the estimators have lower rank and enjoy better prediction performance) over nuclear norm based estimates, on a variety of problem instances.

Some recent works (Jain et al., 2010, 2013; Hardt, 2014; Hardt and Wootters, 2014; Chen and Wainwright, 2015; Ma et al., 2017; Chen et al., 2019a) study the scope of alternating minimization or (projected) gradient stylized algorithmic strategies for the rank constrained optimization problem, similar to Problem (1) — see also Hastie et al. (2016) for related discussions. We should emphasize that our work herein, studies the entire family of nonconvex spectral penalized problems of the form of Problem (3), and is hence more general than the class of estimation problems considered in those works. We establish empirically that this flexible family of nonconvex penalized estimators leads to solutions with better statistical properties than those available from particular instantiations of the penalty function — nuclear norm regularization (2) and rank regularization (1). Along the lines of the aforementioned works, there exists an active stream of research on characterizing the global optimality of local algorithms for various matrix factorization based formulations (Bhojanapalli et al., 2016; Ge et al., 2016; Sun and Luo, 2016; Zheng and Lafferty, 2016; Ge et al., 2017; Shapiro et al., 2018). Our paper focuses on a more general family of nonconvex regularization, with admittedly less strong algorithmic guarantees. Finally, a series of iterative reweighted algorithms have been proposed and discussed (Mazumder et al., 2010; Mohan and Fazel, 2010; Fornasier et al., 2011; Mohan and Fazel, 2012; Gu et al., 2017), largely motivated by the reweighting ideas from sparse recovery problems (Zou, 2006; Candes et al., 2008; Daubechies et al., 2010). Different weight formulas have been suggested to improve the statistical and computational efficiency. These are, however, beyond the scope of the current paper.

1.1 Contributions and Outline

The main contributions of our paper can be summarized as follows:

• We propose a computational framework for nonconvex penalized matrix completion problems of the form (3). Our algorithm: NC-Impute, may be thought of as a novel adaptation (with important enhancements and modifications) of the EM-stylized procedure Soft-Impute

• We present an in-depth investigation of nonconvex spectral thresholding operators, which form the main building block of our algorithm. We also study their effective degrees of freedom (df ),

which provide a simple and intuitive way to calibrate the two-dimensional grid of tuning parameters, extending the scope of the method proposed in nonconvex penalized (least squares) regression by Mazumder et al. (2011) to spectral thresholding operators. We propose computationally efficient methods to approximate the df using tools from random matrix theory.

• We provide comprehensive computational guarantees of our algorithm — this includes the number of iterations needed to reach a first order stationary point and the asymptotic convergence of the sequence of estimates produced by NC-Impute.

• Every iteration of NC-Impute requires the computation of a low-rank SVD of a structured matrix, for which we propose new methods. Using efficient warm-start tricks to speed up the low-rank computations, we demonstrate the effectiveness of our proposal to large scale instances up to the Netflix size in reasonable computation times.

• Over a wide range of synthetic and real-data examples, we show that our proposed nonconvex penalized framework leads to high quality solutions with excellent statistical properties, which are often found to be significantly better than nuclear norm regularized solutions in terms of producing low-rank solutions with good predictive performances.

• Implementations of our algorithms in the R programming language have been made publicly available on github at: https://github.com/diegofrasal/ncImpute.

The remainder of the paper is organized as follows. Section 2 studies several properties of nonconvex spectral penalties and associated spectral thresholding operators, including their effective degrees of freedom. Section 3 describes our algorithmic framework NC-Impute and studies the convergence properties of the algorithm. Section 4 presents numerical experiments demonstrating the usefulness of nonconvex penalized estimation procedures in terms of superior statistical properties on several synthetic datasets — we also show the usefulness of these estimators on several real data instances. Section 5 contains the conclusions and discusses several important future research directions. To improve readability, some technical materials and empirical results are relegated to Section 6.

Notation: For a matrix , we denote its (i, j)th entry by ) is a matrix with its (i, j)th entry given by Ω and zero otherwise, with Ω use the notation ) to denote the projection of A onto the complement of Ω. Let denote the singular values of ) – we will use the notation ) to denote the vector of singular values. When clear from the context, we will simply write instead of ). For a vector , we will use the notation diag(a) to denote an diagonal matrix with ith diagonal entry being

We begin our analysis by considering the fully observed version of Problem (3), given by:

where, for a given matrix Z, a minimizer of the function g(X), denoted by thresholding operator induced by the spectral penalty the SVD of Z. For the nuclear norm regularized problem with the penalty function the corresponding thresholding operator, denoted by ) (say), is given by the familiar soft-thresholding operator (Cai et al., 2010; Mazumder et al., 2010):

where, th entry of separability of the thresholding operator). Here, ) is the the soft-thresholding operator on the singular values of Z and plays a crucial role in the Soft-Impute algorithm (Mazumder et al., 2010). For the rank regularized problem, with

the thresholding operator denoted by ) is given by the hard-thresholding operator (Mazumder et al., 2010):

with A closely related thresholding operator that retains the top r singular values and sets the remaining to zero formed the basis of the Hard-Impute algorithm in Mazumder et al. (2010); Troyanskaya et al. (2001). The results in (5) and (6) suggest a curious link — the spectral thresholding operators (for the two specific choices of the spectral penalty functions given above) are tied to the corresponding thresholding functions that operate only on the singular values of the matrix — in other words, the operators change the singular vectors of the matrix Z. It turns out that a similar result holds true for more general spectral penalty functions as the following proposition illustrates.

denote the SVD of denote the following thresholding operator on the singular values of Z:

Proof. Note that by the Wielandt-Hoffman inequality (Horn and Johnson, 2012) we have that: where, for a vector denotes the standard Euclidean norm. Equality holds when X and Z share the same left and right singular vectors. This leads to:

In the above inequality, note that the left hand side is g(X) (defined in (4)) and right hand side is ¯)) (defined in (7)). It follows that

where, we used the observation that as defined in (7) minimizes ¯addition, this minimum is attained by the function g(X), at the choice completes the proof of the proposition.

Due to the separability of the optimization Problem (7) across the coordinates, i.e., ¯) (where, ¯) is defined in (9)), it suffices to consider each of the subproblems separately. Let ) denote a minimizer of ¯

It is easy to see that the ith coordinate of ) is given by ). This discussion suggests that our understanding of the spectral thresholding operator ) is intimately tied to the univariate thresholding operator (9). Thusly motivated, in the following, we present a concise discussion about univariate penalty functions and the resultant thresholding operators. We begin with some examples of concave penalties that are popularly used in statistics in the context of sparse linear modeling.

Families of Nonconvex Penalty Functions: Several types of nonconvex penalties are popularly used in high-dimensional regression frameworks—see for example, Nikolova (2000); Lv and Fan (2009); Zhang and Zhang (2012). For our setup, since these penalty functions operate on the singular values of a matrix, it suffices to consider nonconvex functions that are defined only on the nonnegative real numbers. We present a few examples below:

• The penalty (Frank and Friedman, 1993) given by

• The SCAD penalty (Fan and Li, 2001) is defined via:

where ) denotes the derivative of

• The MC+ penalty (Zhang, 2010; Mazumder et al., 2011) defined as

Figure 1 shows some members of the above nonconvex penalty families. The penalty function is non differentiable at = 0, due to the unboundedness of 0+. The nonzero derivative at = 0+ encourages sparsity. The penalty functions show a clear transition from the penalty to the penalty — similarly, the resultant thresholding operators show a passage from the soft-thresholding to the hard-thresholding operator. Let us examine the analytic form of the thresholding function induced by the MC+ penalty (for any

Figure 1: [Top panel] Examples of nonconvex penalties = 1 for different values of . [Bottom Panel] The corresponding scalar thresholding operators: = 1, some of the thresholding operators corresponding to the penalty function are discontinuous, and some of the other thresholding functions are “close” to being so.

It is interesting to note that for the MC+ penalty, the derivatives are all bounded and the thresholding functions are continuous for all , the threshold operator (10) coincides with the soft-thresholding operator. However, as 1+ the threshold operator approaches the discontinuous hard-thresholding operator ) — this is illustrated in Figure 1 and can also be observed by inspecting (10). Note that the penalty penalizes small and large singular values in a similar fashion, thereby incurring an increased bias in estimating the larger coefficients. For the MC+ and SCAD penalties, we observe that they penalize the larger coefficients less severely than the simultaneously, they penalize the smaller coefficients in a manner similar to that of the On the other hand, the penalty (for small values of ) imposes a more severe penalty for values of 0, quite different from the behavior of other penalty functions. In general, for a given family of nonconvex penalties ), the effect of () on the nonconvexity can be characterized through the general concavity quantity that is to be introduced in (11).

2.1 Properties of Spectral Thresholding Operators

The nonconvex penalty functions described in the previous section are concave functions on the nonnegative real line. We will now discuss measures that may be thought (loosely speaking) to measure the amount of concavity in the functions. For a univariate penalty function assumed to be differentiable on (0), we introduce the following quantity () that measures the amount of concavity (see also, Zhang (2010)) of

where ) denotes the derivative of

holds:

for some . In inequality (12), subgradient (assuming it exists) of = 0 then the function is simply convex. Using standard properties of spectral functions (Borwein and Lewis, 2006; Lewis, 1995), it follows that -strongly convex iff the vector function:

is -strongly convex on ) denotes the singular values of Z. Let us recall the separable decomposition of ¯) as defined in (9). Clearly, the function -strongly convex (on the nonnegative reals) iff each summand ¯-strongly convex on 0. Towards this end, notice that ¯) is convex on 0 — in particular, ¯is -strongly convex with parameter , provided this number is nonnegative. In this vein, we have the following proposition:

, then the function -strongly convex with

For the MC+ penalty, the condition 0 is equivalent to 1. For the function, with 1, the parameter , and thus the function g(X) is not strongly convex.

is Lipschitz continuous with constant i.e, for all

Proof. We rewrite g(X) as:

), and rearranging the terms in (15), it follows that ), a minimizer of g(X), is given by:

If 0, the function ) is convex for every 0, then the first term appearing in the objective function in (16) is convex. Thus, assuming 0 both summands in the above objective function are convex. In particular, the optimization problem (16) is convex and ) can be viewed as a convex proximal map (Rockafellar, 1970). Using standard contraction properties of proximal maps, we have that:

Since the above holds true for any as chosen above, optimizing over the value of Problem (16) remains convex gives us ˆ), thereby leading to (14).

2.2 Effective Degrees of Freedom for Spectral Thresholding Operators

In this section, to better understand the statistical properties of spectral thresholding operators, we study their degrees of freedom. The effective degrees of freedom or df is a popularly used statistical notion that measures the amount of “fitting” performed by an estimator (Efron et al., 2004; Hastie et al., 2009; Stein, 1981). In the case of classical linear regression, for example, df is simply given by the number of features used in the linear model. This notion applies more generally to additive fitting procedures. Following Efron et al. (2004); Stein (1981), let us consider an additive model of the form:

for ), for the fully observed model above, i.e., (17) is given by:

where denotes the (i, j)th entry of the matrix . For the particular case of a spectral thresholding operator we have ˆ) satisfies a weak differentiability condition, the df may be computed via a divergence formula (Stein, 1981; Efron et al., 2004):

where (For the spectral thresholding operator ), expression (18) holds if the map ) is Lipschitz and hence weakly differentiable — see for example, Cand`es et al. (2013). In the light of Proposition 3, the map ) is Lipschitz when Under the model (17), the singular values of Z will have a multiplicity of one with probability one. We assume that the univariate thresholding operators are differentiable, i.e., ) exists. With these assumptions in place, the divergence formula for ) can be obtained following Cand`es et al. (2013); Mazumder and Weng (2020), as presented in the following proposition.

Proposition 4. Assume that 1 + and the model (17) is in place. Then the degrees of freedom of the estimator is given by:

where the ’s are the singular values of Z.

We note that the above expression is true for any value of 1 + 0. For the MC+ penalty function, expression (19) holds for 1. As soon as 1, the above method of deriving df does not apply due to the discontinuity in the map ). Values of close to one (but larger),

Figure 2: Figure showing the df for the MC+ thresholding operator for a matrix with and v = 1. The df profile as a function of (in the log scale) is shown for three values of . The dashed lines correspond to the df of the spectral soft-thresholding operator, corresponding to . We propose calibrating the () grid to a () grid such that the df corresponding to every value of matches the df of the soft-thresholding operator — as shown in Figure 3.

however, give an expression for the df near the hard-thresholding spectral operator, which corresponds to

To understand the behavior of the df as a function of (), let us consider the null model with = 0 and the MC+ penalty function. In this case, for a fixed (see Figure 2 with a fixed df is seen to increase with smaller values: the soft-thresholding function shrinks the large coefficients and sets all coefficients smaller than to be zero; the more aggressive (closer to the hard thresholding operator) shrinkage operators ()) shrink less for larger values of and set all coefficients smaller than to zero. Thus, intuitively, the more aggressive thresholding operators should have larger df since they do more “fitting” — this is indeed observed in Figure 2. Mazumder et al. (2011) studied the df of the univariate thresholding operators in the linear regression problem, and observed a similar pattern in the behavior of the ) values. For the linear regression problem, Mazumder et al. (2011) argued that it is desirable to choose a parametrization for () such that for a fixed , as one moves across should be the same. We follow the same strategy for the spectral regularization problem considered herein — we reparametrize a two-dimensional grid of (to a two-dimensional grid of () values, such that the df remain calibrated in the sense described above — this is illustrated in Figure 2 (see the horizontal dashed lines corresponding to the constant df values, after calibration). Figure 3 shows the lattice of () after calibration. The values of (˜on each curve induce the same moves down (the penalty becomes more “nonconvex”), the corresponding ˜(the shrinkage) has to increase to maintain the same df.

The study of df presented herein provides a simple and intuitive explanation about the roles played by the parameters () for the fully observed problem. The notion of calibration provides a new parametrization of the family of penalties. From a computational viewpoint, since the general algorithmic framework presented in this paper (see Section 3) computes a regularization surface using warm-starts across adjacent () values on a two-dimensional grid; it is desirable for the adjacent values to be close — the df calibration also ensures this in a simple and intuitive manner.

Computation of df : The df estimate as implied by Proposition 4 depends only upon singular values (and not the singular vectors) of a matrix and can hence be computed with cost The expectation can be approximated via Monte-Carlo simulation — these computations are easy to

Figure 3: Figure showing the calibrated () lattice — for every fixed value of of the MC+ spectral threshold operators are the same across different values. The df computations have been performed on a null model using Proposition 5.

parallelize and can be done offline. Since we compute the df for the null model, for larger values of m, n we recommend using the Marchenko-Pastur law for iid Gaussian matrix ensembles to approximate the df expression (19). We illustrate the method using the MC+ penalty for 1. Towards this end, let us define a function on

For the following proposition, we will assume (for simplicity) that

, then under the model

where is the thresholding operator corresponding to the MC+ penalty with the expectation is taken with respect to independently generated from the Marchenko-Pastur distribution (see Lemma 1, Section 6.1).

Proof. For a proof, see Section 6.1.1.

Note that the variance vin model (17) can always be assumed to be one (by adjusting the value of the tuning parameter accordingly

In this section, we present algorithm NC-Impute. The algorithm is inspired by an EM-stylized procedure, similar to Soft-Impute (Mazumder et al., 2010), but has important innovations, as we will discuss shortly. It is helpful to recall that, for observed data: ), the algorithm SoftImpute relies on the following update sequence

which can be interpreted as computing the nuclear norm regularized spectral thresholding operator for the following “fully observed” problem:

where, the missing entries are filled in by the current estimate, i.e., ). We refer the reader to Mazumder et al. (2010) for a detailed study of the algorithm. Mazumder et al. (2010) suggest, in passing, the notion of extending Soft-Impute to more general thresholding operators; however, such generalizations were not pursued by the authors. In this paper, we present a thorough investigation about nonconvex generalized thresholding operators — we study their convergence properties, scalability aspects and demonstrate their superior statistical performance across a wide range of numerical experiments. Update (20) suggests a natural generalization to more general nonconvex penalty functions, by simply replacing the spectral soft thresholding operator ) with more general spectral operators

While the above update rule works quite well in our numerical experiments, it enjoys limited computational guarantees, as suggested by our convergence analysis in Section 3.1. We thus propose and study a seemingly minor generalization of the rule (21) — this modified rule enjoys superior finite time convergence rates to a first order stationary point. We develop our algorithmic framework below.

Let us define the following function:

for 0. Note that ) majorizes the objective function f(X) defined in (3), i.e., , with equality holding at . In an attempt to obtain a minimum of Problem (3), we propose to iteratively minimize ), an upper bound to f(X), to obtain — more formally, this leads to the following update sequence:

Note that is easy to compute; by some rearrangement of (22) we see:

where + 1). Note that (24) is a minor modification of (21) — in particular, if = 0, then these two update rules coincide.

The sequence defined via (24) has desirable convergence properties, as we discuss in Section 3.1. In particular, as , the sequence reaches (in a sense that will be made more precise later) a first order stationary point for Problem (3). We also provide a finite time convergence analysis of the update sequence (24).

We intend to compute an entire regularization surface of solutions to Problem (3) over a two-dimensional grid of ()-values, using warm-starts. We take the MC+ family of functions as a running example, with (beginning, we compute a path of solutions for the nuclear norm penalized problem, i.e., Problem (3) with on a grid of values. For a fixed value of , we compute solutions to Problem (3) for smaller values of , gradually moving away from the convex problems. In this continuation scheme, we found the following strategies useful:

• For every value of (), we apply two copies of the iterative scheme (23) initialized with solutions obtained from its two neighboring points (). From these two candidates, we select the one that leads to a smaller value of the objective function (

• Instead of using a two-dimensional rectangular lattice, one can also use the recalibrated lattice, suggested in Section 2.2, as the two-dimensional grid of tuning parameters.

The algorithm outlined above, called NC-Impute is summarized as Algorithm 1.

We now present an elementary convergence analysis of the update sequence (24). Since the problems under investigation herein are nonconvex, our analysis requires new ideas and techniques beyond those used in Mazumder et al. (2010) for the convex nuclear norm regularized problem.

3.1 Convergence Analysis

By the definition of we have that:

Let us define the quantities:

where, if 0, then the function )-strongly convex. In particular, from (23), it follows that ), a subgradient of the map ) (evaluated at zero. We thus have:

Now note that, by the definition of , we always have: which combined with (25) leads to (replacing

In addition, we have:

Combining (26) and (27), and observing that ), we have:

Since ∆0, the above inequality immediately implies that and the improvement in objective values is at least as large as the quantity ∆). The term ∆) is a measure of progress of the algorithm, as formalized by the following proposition.

and for any , let us consider the update Then the following are equivalent:

(iii) is a fixed point, i.e.,

= 0, we have that . The condition implies that

where the term on the right equals is a fixed point.

say — this implies that ∆. Let us now consider two cases, depending upon the value of 0, then we have . On the other hand, if the quantities = 0, the conclusion needs to be modified: ∆0 implies that

Motivated by the above discussion, we make the following definition of a first order stationary point for Problem (3).

is said to be a first order stationary point for Problem (3) if ∆is said to be an -accurate first order stationary point for Problem (3) if ∆

Proposition 7. The sequence is decreasing and suppose it converges to ˆf. Then the rate of convergence of to this first order stationary point is given by:

Proof. The arguments presented preceding Proposition 7 establish that the sequence creasing and converges to ˆf, say. Consider (28) for any 1 . We have that ∆) — summing this inequality for k = 1, . . . , K we obtain:

where in the last inequality we used the simple fact that . Gathering the left and right parts of the above chain of inequalities leads to (29).

Proposition 7 shows that the sequence reaches an -accurate first order stationary point within many iterations. The number of iterations , depends upon how close the initial estimate ) is to the eventual solution ˆf. Since NC-Impute employs warm-starts, the constant appearing in the rhs of (29) suggests that the number of iterations required to a reach an approximate first order stationary point is quite low — this is indeed observed in our experiments, and this feature of using warm-starts makes our algorithm particularly attractive from a practical viewpoint.

3.1.1 Rank Stabilization

Let us consider the thresholding function ) defined in (24), which expresses as a function of . Using the development in Section 2, it is easy to see that the spectral operator closely tied to the following vector thresholding operator (30), acting on the singular values of Formally, for a given nonnegative vector , if we denote:

where is the SVD of Thus, properties of the thresholding function are closely related to those of the vector thresholding operator ). Due to the separability of the vector thresholding operator ), across each coordinate of , we denote by coordinate of

We now investigate what happens to the rank of the sequence as defined via (23). In particular, does this rank converge? We show that the rank stabilizes after finitely many iterations, under an additional assumption — namely the spectral thresholding operator is discontinuous — see Figure 1 for examples of discontinuous thresholding functions.

Proposition 8. Consider the update sequence

as defined in (24); and let . Suppose that there is a such that, for any scalar the following holds: — i.e., the scalar thresholding operator is discontinuous. Then there exists an integer such that for all i.e., the rank stabilizes after finitely many iterations.

Proof. Using (28) it follows that

where the last inequality follows from Wielandt-Hoffman inequality (Horn and Johnson, 2012) and ) denotes the vector of singular values of ) be an indicator vector with ith coordinate being equal to 1(= 0). We will prove the result of rank stabilization via the method of contradiction. Suppose the rank does not stabilize, then ) for infinitely many k values. Thus there are infinitely many values such that:

where i is taken such that = 0. Note that by the property of the thresholding function ) we have that . This implies that infinitely many values, which is a contradiction to the convergence: 0. Thus the support of ) converges, and necessarily after finitely many iterations — leading to the existence of an iteration number , after which the rank of remains fixed. This completes the proof of the proposition.

, the discontinuity of the thresholding operator (as demanded by Proposition 8) occurs for the MC+ penalty function as soon as . For a general , discontinuity in occurs as soon as

3.1.2 Subspace Stabilization

We study herein, the properties of the left and right singular subspaces associated with the sequence of subspaces has important implications in the main bottleneck of the NCImpute algorithm, i.e., the SVD computations — we discuss this in further detail in Section 3.2. The study of singular subspace stabilization requires subtle analysis based on matrix perturbation theory (Stewart and Sun, 1990), since the left (and right) singular subspace, corresponding to the top r singular values of a matrix is not a continuous function of the matrix argument.

Towards this end, we first recall a standard notion of distance between two subspaces (with same dimension) in terms of canonical angles.

be two orthonormal matrices and let us define that [forms an orthonormal basis for . The canonical angles between these two subspaces denoted by the vector Θ(are defined as:

where, are the singular values of the matrix

We now present a result regarding perturbation of singular subspaces, taken from Stewart and Sun (1990). Before stating the proposition, we introduce some notation. Let denote a matrix of the left singular vectors (respectively, right) of a matrix diagonal matrix of the corresponding top singular values. Similarly, we use the notation to denote the triplet of left and right singular vectors and singular values (corresponding to the top singular values) for a matrix . We use the following matrices

to measure a notion of proximity between . The distance between the left (and also right) singular subspaces (corresponding to the top singular values) of may be measured by the following quantity:

where, the notation denotes the spectral norm of A. With the above notations in place, we present the following proposition (Stewart and Sun, 1990) regarding the perturbation of singular subspaces of matrices.

Proposition 9. Suppose there exists

where, Σis a diagonal matrix with the remaining singular values of A. Then,

The above proposition informs us about the proximity of the left (and also right) singular subspaces across successive iterates , as presented in the following proposition:

Proof. The proof is presented in Section 6.1.2.

Remark 2. Let the assumptions of Proposition 8 be in place – this implies that there exists an integer

Hence, in particular, there is a separation between sufficiently large. This implies that , i.e., in words: the distance between the left (and right) singular subspaces corresponding to the top r singular values of converges to zero, as

3.1.3 Asymptotic Convergence.

We now investigate the asymptotic convergence properties of the sequence Proposition 8 shows that under suitable assumptions, the sequence rank(1 converges. The existence of a limit point of is guaranteed if the singular values of ) remain bounded. It is not immediately clear whether the sequence ) will remain bounded since several spectral penalty functions (like the MC+ penalty) are bounded. We address herein, the existence of a limit point of the sequence ), and hence the sequence

For the following proposition, we will assume that the concave penalty function 0 is differentiable and the gradient is bounded.

denote the rank-reduced SVD of a limit point of the sequence , such that (along a subsequence denote the ith column of ¯U (and similarly for ¯) and let us denote ¯Θ = [vec(. We have the following:

(b) If , then the sequence converges to a first order stationary point: ¯

Proof. See Section 6.1.3

Proposition 8 describes sufficient conditions under which the rank of the sequence stabilizes after finitely many iterations — it does not describe the boundedness of the sequence is addressed in Proposition 11. Note that Proposition 11 does not imply that the rank of the sequence stabilizes after finitely many iterations (recall that Proposition 11 does not assume that the thresholding operators are discontinuous, an assumption required by Proposition 8).

3.2 Computing the Thresholding Operators

The operator (24) requires computing a thresholded SVD of the matrix , as demonstrated by Proposition 1. The thresholded singular values ) as in (30) will have many zero coordinates due to the “sparsity promoting” nature of the concave penalty. Thus, computing the thresholding operator (24) will typically require performing a low-rank SVD on the matrix While direct factorization based SVD methods can be used for smaller problems where min{m, n} is of the order of a thousand or so; for larger matrices, such methods become computationally prohibitive — we thus resort to iterative methods for computing low-rank SVDs for large scale problems. Algorithms such as the block power method; also known as block QR iterations, or those based on the Lanczos method (Golub and Van Loan, 1983) are quite effective in computing the top few singular value and vectors of a matrix A, especially when the operations of multiplying (for vectors of matching dimensions) can be done efficiently. Indeed, such matrix-vector multiplications turn out to be quite computationally attractive for our problem, since the computational cost of multiplying

which admits a decomposition as the sum of a sparse matrix and a low-rank matrix. Note that the sparse matrix has the same sparsity pattern as the observed indices Ω. Decomposition (32) is inspired by a similar decomposition that was exploited effectively in the algorithm Soft-Impute(Mazumder et al., 2010), where the authors use PROPACK (Larsen, 2004) to compute the low-rank SVDs. In this paper, we use the Alternating Least Squares (ALS)-stylized procedure, which computes a low-rank SVD by solving the following nonlinear optimization problem:

using alternating least squares—this is in fact, equivalent to the block power method (Golub and Van Loan, 1983), in computing a rank ˜r SVD of the matrix . Across the iterations of NC-Impute, we pass the warm-start information in the U, V ’s obtained from a low-rank SVD of to compute the low-rank SVD for . Empirically, this warm-start strategy is found to be significantly more advantageous than a black-box low-rank SVD stylized approach, as used in the Soft-Impute algorithm (for example), where, at every iteration, a new low-rank SVD is computed from scratch via PROPACK. This strategy quite naturally leads to a loss of useful information about the left and right singular vectors, which become closer to each other along the course of the Soft-Impute iterations (as formalized by Section 3.1.2). Using warm-start information across successive iterations (i.e., k values) leads to notable gains in computational speed (often reduces the total time to compute a family of solutions by orders of magnitude), when compared to black-box SVD stylized methods that do not rely on such warm-start strategies. This improvement is also supported by theory — the computational guarantee of block power iterations (Golub and Van Loan, 1983) states that the subspace spanned by the U matrix (in the factorization in (33)) converges to that of the top ˜r left singular vectors at the rate: denotes the number of power iterations, depends upon the ratio between the ˜singular values of the matrix depends upon the distance between: the initial estimate of (the subspace spanned by) U and the left top-˜r set of singular vectors of ˜. The constant C is smaller with a good warm-start, when compared to a random initialization. A similar argument applies for the right set of singular vectors.

In this section, we present a systematic experimental study of the statistical properties of estimators obtained from (3) for different choices of penalty functions. We perform our experiments on a wide array of synthetic and real data instances. Recall that the majority of the algorithmic guarantees proved in Section 3 rely on the condition 0. For MC+ penalty functions, it is straightforward to verify that 0 as long as ]. Hence we will use throughout this section.

4.1 Synthetic Examples

We study three different examples, where, for the true low-rank matrix , we vary both the structure of the left and right singular vectors in L and R, as well as the sampling scheme used to obtain the observed entries in Ω. Our basic model is , where we observe entries

Figure 4: (Color online) Random Orthogonal Model (ROM) simulations with SNR = 1. The choice refers to nuclear norm regularization as provided by the Soft-Impute algorithm. We also include the choice = 1 to represent the rank regularized approach. The least nonconvex alternatives at behave similarly to nuclear norm, although with better prediction performance. The choices of in excessively aggressive fitting behavior for the true rank = 10 case, but improve significantly in prediction error and recovering the true rank in the sparser true rank = 5 setting. In both scenarios, the intermediate models with = 20 fare the best, with the former achieving the smallest prediction error, while the latter estimates the actual rank of the matrix. Values of test error larger than one are not displayed in the figure.

(Ω. We consider different types of missing patterns for Ω, and various signal-to-noise (SNR) ratios for the Gaussian error term , defined here to be:

Accordingly, the (standardized) training and test error for the model are defined as:

where a value greater than one for the test error indicates that the computed estimate ˆM does a worse job at estimating M than the zero solution, and the training error corresponds to the fraction of the error explained on the observed entries by the estimate ˆM relative to the zero solution.

Figure 5: (Color online) Random Orthogonal Model (ROM) simulations with SNR = 5. The benefits of nonconvex regularization are more evident in this high-sparsity, high-missingness scenario. While the and = 80 models distance themselves more from nuclear norm, the remaining members of the MC+ family essentially minimize prediction error while correctly estimating the true rank. This is especially true in panel (d), where the best predictive performance of the model = 5 at the correct rank is achieved under a low-rank truth and high SNR setting.

Example-A: In our first simulation setting, we use the model

where L and R are matrices generated from the random orthogonal model (Cand`es and Recht, 2009), and the singular values Φ = diag() are randomly selected as The set Ω is sampled uniformly at random. Recall that for this model, exact matrix completion in the noiseless setting is guaranteed as long as , for some universal constant C (Recht, 2011). Under the noisy setting, Mazumder et al. (2010) show superior performance of nuclear norm regularization vis-`a-vis other matrix recovery algorithms (Cai et al., 2010; Keshavan et al., 2010) in terms of achieving smaller test error. For the purposes herein, we fix (m, n) = (800, 400) and set the fraction of missing entries to

Example-B: In our second setting, we also consider the model

Figure 6: (Color online) Coherent and Nonuniform Sampling (NUS) simulations with SNR = 10. nonconvex regularization also proves to be a successful strategy in these challenging scenarios, particularly in the nonuniform sampling setting where the MC+ family exhibits a monotone decrease in prediction error as approaches 1. Again, the model = 5 estimates the correct rank under high SNR settings. Although nuclear norm achieves a relatively small prediction error, compared with previous simulation settings, the MC+ family still provides a superior and more robust mechanism for regularization.

but we now select matrices L and R which do not satisfy the incoherence conditions required for full matrix recovery. Specifically, for the choices of (9, we select L and R to be block-diagonal matrices of the form

where 5, are random matrices with scaled Gaussian entries. The singular values are again sampled as 100) with Ω being uniformly random over the set of indices. For this model, successful matrix completion is not guaranteed even for the noiseless problem with the nuclear norm relaxation, as the left and right singular vectors are not sufficiently spread. We observe the usefulness of the nonconvex regularized estimators in this regime, in our experimental results.

Example-C: In our third simulation setting, for the choice of (m, n, r) = (100, 100, 10), we also generate from the random orthogonal model as in our first setting, but we now allow the observed entries in Ω to follow a nonuniform sampling scheme. In particular,

Figure 7: (Color online) MovieLens 100k and 1m data. For each value of in the solution path, an operating rank threshold (capped at 250) larger than the rank of the previous solution was employed.

we fix Ω

with the fraction of missing entries thus being 25. This is again a challenging simulation setting in which both the uniform (Cand`es and Recht, 2009) and independent (Chen et al., 2014) sampling scheme assumptions in Ω are violated. Our aim again is to explore whether the nonconvex MC+ family is able to outperform nuclear norm regularization in this regime.

For all three settings above, we choose a 10025 grid of () values as follows. In each simulation instance we fix , the smallest value of for which the nuclear norm regularized solution is zero, and set . Keeping in mind that NC-Impute benefits greatly from using warm starts, we construct an equally spaced sequence of 100 values of decreasing from choose 25 -values in a logarithmic grid from 5000 to 1.1. The results displayed in Figures 4 – 6 show averages of training and test errors, as well as recovered ranks of the solution matrix ˆvalues of (), taken over 50 simulations under all three problem instances. The plots including rank reveal how effective the MC+ family is at recovering the true rank while minimizing prediction error. Throughout the simulations we keep an upper bound of the operating rank as 50.

4.1.1 Discussion of Experimental Results

We devote Figures 4 and 5 to analyze the simpler random orthogonal model (Example-A), leaving the more challenging coherent and nonuniform sampling settings (Example-B and Example-C) for Figure 6. In each case, the captions detail the results which we summarize here. The noise is quite high in Figure 4 with SNR = 1 and 90% of the entries missing in both displayed settings, while the model complexity decreases from a true rank of 10 to 5. The underlying true ranks remain the same in Figure 5, but the noise level has decreased to SNR = 5 with the missing entries increasing to 95%. For each model setting considered, all nonconvex methods from the MC+ family outperform nuclear norm regularization in terms of prediction performance, while members of the MC+ family with smaller values of are better at estimating the correct rank. The choices of have the best performance in Figure 4 (best prediction errors around the true ranks), while more nonconvex alternatives fare better in the high-sparsity, low-noise setting of Figure 5. In both figures, the performance of nuclear norm regularization is somewhat similar to the least nonconvex alternative displayed at = 100, however, the bias induced in the estimation of the singular values of the low-rank matrix M leads to the worst bias-variance trade-off among all training versus test error plots for the settings considered.

While the nuclear norm relaxation provides a good convex approximation for the rank of a matrix (Recht et al., 2010), these examples show that nonconvex regularization methods provide a superior mechanism for rank estimation. This is reminiscent of the performance of the MC+ penalty in the context of variable selection within high-dimensional sparse regression models. Although the penalty function represents the best convex approximation to the penalty, the gap bridged by the nonconvex MC+ penalty family provides a better basis for model selection, and hence rank estimation in the low-rank matrix completion setting.

For the coherent and nonuniform sampling settings of Figure 6, we choose the small noise scenario SNR = 10 in order to favor all considered models. Despite the absence of any theoretical guarantees for successful matrix recovery, the nuclear norm regularization approach achieves a relatively small prediction error in all displayed instances. Nevertheless, the nonconvex MC+ family of penalties seems empirically more adept at overcoming the limitations of nuclear norm penalizated matrix completion in these challenging simulation settings. In particular, the most aggressive nonconvex fitting behavior at = 5 achieves excellent prediction performance in the nonuniform sampling setting while correctly estimating the true rank of the coherent model.

4.2 Real Data Examples: MovieLens and Netflix datasets

We now use the real world recommendation system datasets ml100k and ml1m provided by MovieLensas well as the famous Netflix competition data to compare the usual nuclear norm approach with the MC+ regularizers. The dataset ml100k consists of 100, 000 movie ratings (1–5) from 943 users on 1, 682 movies, whereas ml1m includes 1, 000, 209 anonymous ratings from 6, 040 users on 3, 952 movies. In both datasets, for all regularization methods considered, a random subset of 80% of the ratings were used for training purposes; the remaining were used as the test set.

We also choose a similar 100 25 grid of () values, but for each value of in the decreasing sequence, we use an “operating rank” threshold somewhat larger than the rank of the previous solution, with the goal of always obtaining solution ranks smaller than the operating threshold. Following the approach of Hastie et al. (2016), we perform row and column centering of the corresponding (incomplete) data matrices as a preprocessing step.

Figure 7 compares the performance of nuclear norm regularization with the MC+ family of penalties on these datasets, in terms of the prediction error (RMSE) obtained from the left out portion of the

Figure 8: (Color online) Netflix competition data. The model = 10 achieves optimal test set RMSE of 0.8276 for a solution rank of 105.

data. While the fitting behavior at = 5 seems to be overly aggressive in these instances, the choice = 10 achieves the best test set RMSE with a minimum solution rank of 20 for the ml100k data. With a higher test RMSE, nuclear norm regularization achieves its minimum with a less parsimonious model of rank 62. Similar results hold for the ml1m data, where the model = 15 achieves near optimal test RMSE at a solution rank of 115, while the best estimation accuracy of Soft-Impute occurs for ranks well over 200.

The Netflix competition data consists of 100, 480, 507 ratings from 480, 189 users on 17, 770 movies. A designated probe set, a subset of 1, 408, 395 of these ratings, was distributed to participants for calibration purposes, leaving 99, 072, 112 for training. We did not consider the probe set as part of this numerical experiment, instead choosing 1, 500, 000 randomly selected entries as test data with the remaining 97, 572, 112 used for training purposes. Similar to the MovieLens data, we select a 20 25 grid of () values which adaptively chooses an operating rank threshold, and also remove row and columns means for prediction purposes.

As shown in Figure 8, the MC+ family again yields better prediction performance under more parsimonious models. On average, and for a convergence tolerance of 0.001 in Algorithm 1, the sequence of twenty models took under 10.5 hours of computing on an Intel E5-2650L cluster with 2.6 GHz processor. We note that our main goal here is to show the feasibility of applying NC-Impute to the MC+ family on a Netflix sized dataset, and further reductions in computation time may be possible with specialized implementations. It seems that using a family of nonconvex penalties leads to models with better statistical properties, when compared to the nuclear norm regularized problem and the rank constrained problem (obtained via Hard-Impute, for example).

In this paper we present a computational study for the noisy matrix completion problem with non-convex spectral penalties — we consider a family of spectral penalties that bridge the convex nuclear norm penalty and the rank penalty, leading to a family of estimators with varying degrees of shrinkage and nonconvexity. We propose NC-Impute— an algorithm that appropriately modifies and enhances the EM-stylized procedure Soft-Impute (Mazumder et al., 2010), to compute a two dimensional family of solutions with specialized warm-start strategies. The main computational bottleneck of our algorithm is a low-rank SVD of a structured matrix, which is performed using a block QR stylized strategy that makes effective use of singular subspace warm-start information across iterations. We discuss computational guarantees of our algorithm, including a finite time complexity analysis to a first order stationary point. We present a systematic study of various statistical and structural properties of spectral thresholding functions, which form a main building block in our algorithm. We demonstrate the impressive gains in statistical properties of our framework on a wide array of synthetic and real datasets. The current work leaves open several important directions for future research.

• Statistical guarantees for the nonconvex method. In addition to the comprehensive algorithmic analysis presented in the paper, it is of great importance to establish statistical theories such as estimation error bounds to shed new lights on the empirical success of the proposed nonconvex method. In particular, given that our algorithm returns stationary points, it would be interesting to obtain statistical errors for these local optima. In fact, such type of results have been derived for regularized M-estimators in some general multivariate analysis settings when the penalty function is separable (Loh and Wainwright, 2015). A notable difficulty in the matrix completion problem is that the penalty is imposed on the singular values and is hence a non-separable function of the matrix, which will require more delicate analyses. Moreover, we should point out that statistical analysis of nonconvex optimization methods for matrix completion has been actively investigated in recent years. See the works surveyed in the last paragraph of Introduction and Chi et al. (2019) for a thorough review. However, the nonconvex methods studied in this line of research are exclusively based on low-rank matrix factorization formulation, and the regularization from these methods is less general than the ones in our method. The nonconvexity of the former methods is largely due to the matrix factorization which enables the reduction of memory and computation costs, while the nonconvexity of our method arises from the nonconvex penalties that aim to attenuate the bias.

• Sharp comparison between nonconvex and convex methods. Referring to both simulation study and real data analysis in Section 4, we observe that the value of leading to the optimal matrix completion performance lies between 5 and 30. Recall that the penalty parameter in the MC+ penalties controls the amount of nonconvexity in the regularization. As decreases from to 1, the penalty behaves closer to and farther away from . Hence, the empirical results in Section 4 demonstrate that neither the convex approach () nor the most aggressive nonconvex one (= 1) is the optimal choice. This phenomenon is in fact a manifestation of bias-variance tradeoff. A smaller value of brings more “nonconvexity” to the regularization and hence induces less bias as expected. On the other hand, more “nonconvexity” means more aggressiveness in the selection of low rank matrices and thus results in larger variance. Consequently, for a given level of signal-to-noise ratio in the observations, the optimal is the one that strikes the best balance between bias and variance. This point of view lends further support to our proposed nonconvex method which incorporates the entire family of nonconvex penalties instead of some particular instantiations. Recent works by a subset of the authors (Zheng et al., 2017; Weng et al., 2018; Wang et al., 2019) have given sharp theoretical characterizations of such a phenomenon in the high-dimensional sparse regression and variable selection problems. The results there reveal that among the penalties for 2], as the signal-to-noise ratio (SNR) decreases, the optimal value of q will move from 0 towards 2. See also the work of Hazimeh and Mazumder (2019); Mazumder et al. (2017) for similar observations regarding the overfitting of -based estimators for low SNR regimes. For the MC+ penalties in the matrix completion problem, plays a similar role as q does in the regression problem. It would be of great interest to derive analogue theories for the matrix completion problem and establish a sharp

6.1 Additional Technical Material

Lemma 1. (Marchenko-Pastur law (Bai and Silverstein, 2010)). Let with be the eigenvalues of Define the random spectral measure

Then, assuming

where is a deterministic measure with density

6.1.1 Proof of Proposition 5.

Proof. In the following proof, we make use of the notation: Θ), defined as follows. For two positive sequences ) if there exists a constant c > 0 such that and we say ), whenever,

to represent

Adopting the notation from Lemma 1, it is not hard to verify that

where . A quick check of the relation between

Due to the Lipschitz continuity of the functions ), we obtain

Hence, there exists a positive constant , such that for sufficiently large n,

Let be two independent random variables generated from the Marchenko-Pastur distribution If we can show

then by the Dominated Convergence Theorem (DCT), we conclude the proof in the

regime. Note immediately that

Moreover, given that ) is bounded and continuous, the Marchenko-Pastur theorem in Lemma 1 implies

Since (), and the discontinuity set of the function probability under the measure induced by (), by the continuous mapping theorem,

Also, due to the boundedness of ), it holds that

Combining (34) – (36) completes the proof for the

and

6.1.2 Proof of Proposition 10

Proof. Observe that R as defined in Proposition 9 can be written as:

where, above we have used the fact that , which follows from the definition of the SVD of . By a simple inequality it follows that

where we have used the the fact that = 1. Similarly, we have an analogous result for Q:

6.1.3 Proof of Proposition 11

We set the subdifferential of the map ) to zero at

where is the SVD of . Note that the term: is a subdifferential (Lewis, 1995) of the spectral function:

where is a diagonal matrix with the ith diagonal entry being a derivative of the map 0), denoted by . Note that (40) can be rewritten as:

As , term (a) converges to zero (See Proposition 7) and thus, we have:

Let us denote the ith column of , and use a similar notation for the rank of . Hence, we have:

Multiplying the left and right hand sides of the above by , we have the following:

for denote a limit point of the sequence (which exists since the sequence is bounded); and let r be the rank of ¯. Let us now study the following equations

Using the notation ¯)), we note that (41) are the first order stationary conditions for a point ¯for the following penalized regression problem:

with . If the matrix ¯Θ = [¯] (note that ¯Θ satisfies (41) is finite — in particular, the sequence is bounded and has a limit point: ¯

satisfies the first order stationary condition (41).

then (42) admits a unique solution ¯, which implies that has a unique limit point, and hence the sequence necessarily converges.

6.2 Additional Simulation Results

This section contains additional numerical results from the simulation study in Section 4.1.

• To demonstrate the variation of the procedures in the experiments, we plot the averaged value and standard error of both test error and rank for some representative nonconvex penalty functions. Specifically, under each scenario considered in Section 4.1, we pick the nonconvex penalty that yields the best prediction and rank estimation performance. For each picked penalty, we plot the averaged value of test error and rank along with the associated standard error, against the tuning parameter . The results are shown in Figures 9, 10, and 11. As is clear form the figures, the standard error is typically (at least) one order of magnitude smaller than the average. Moreover, the general patterns of test error and rank on the solution path are expected, except for a few points corresponding to very small values of . The irregularity of these few points occurs probably because the solutions are getting unstable as the nonconvex regularization becomes weak when is significantly small.

• To examine the rank dynamics of the updates in NC-Impute, we compute the number of iterations that the algorithm takes for the convergence of the rank. We choose the same six non-convex penalties as above and evaluate the rank stabilization for several values of . The results are summarized in Figure 12. One clearly observes that except for few instances, it takes less than 10 iterations for the rank to stabilize. Moreover, when the penalty is more “nonconvex” (i.e., is smaller), the rank stabilization occurs earlier. These empirical results provide complementary information on rank stabilization that has been theoretically investigated in 3.1.1.

P. Alquier. A bayesian approach for noisy matrix completion: Optimal rate under general sampling distribution. Electronic Journal of Statistics, 9(1):823–841, 2015.

Z. Bai and J. W. Silverstein. Spectral Analysis of Large Dimensional Random Matrices. Springer, 2010.

D. Bertsimas, A. King, and R. Mazumder. Best subset selection via a modern optimization lens. Annals of Statistics (to appear), 2016.

S. Bhojanapalli, B. Neyshabur, and N. Srebro. Global optimality of local search for low rank matrix recovery. In Advances in Neural Information Processing Systems, pages 3873–3881, 2016.

J. Borwein and A. Lewis. Convex Analysis and Nonlinear Optimization. Springer, 2006.

Figure 9: Random Orthogonal Model (ROM) simulations with SNR = 1. The optimal nonconvex penalties are obtained at = 20 under the two scenarios respectively. The integers from 1 to 100 on the x-axis index the grid of 100 values of (from largest to smallest) as described in Section 4.1.

J.-F. Cai, E. J. Cand`es, and Z. Shen. A singular value thresholding algorithm for matrix completion. SIAM Journal on Optimization, 20:1956–1982, 2010.

E. Cand`es, C. Sing-Long, and J. D. Trzasko. Unbiased risk estimates for singular value thresholding and spectral estimators. Signal Processing, IEEE Transactions on, 61(19):4643–4657, 2013.

E. J. Cand`es and Y. Plan. Matrix completion with noise. Proceedings of the IEEE, 98:925–936, 2010.

E. J. Cand`es and B. Recht. Exact matrix completion via convex optimization. Foundations of Computational mathematics, 9:717–772, 2009.

E. J. Cand`es and T. Tao. The power of convex relaxation: near-optimal matrix completion. IEEE Transactions on Information Theory, 56:2053–2080, 2010.

E. J. Candes, M. B. Wakin, and S. P. Boyd. Enhancing sparsity by reweighted minimization. Journal of Fourier analysis and applications, 14(5-6):877–905, 2008.

Figure 10: Random Orthogonal Model (ROM) simulations with SNR = 5. The optimal nonconvex penalties are obtained at = 5 under the two scenarios respectively. The integers from 1 to 100 on the x-axis index the grid of 100 values of (from largest to smallest) as described in Section 4.1.

J. Chen, D. Liu, and X. Li. Nonconvex rectangular matrix completion via gradient descent without regularization. arXiv preprint arXiv:1901.06116, 2019a.

Y. Chen. Incoherence-optimal matrix completion. IEEE Transactions on Information Theory, 61(5): 2909–2923, 2015.

Y. Chen and M. J. Wainwright. Fast low-rank estimation by projected gradient descent: General statistical and algorithmic guarantees. arXiv preprint arXiv:1509.03025, 2015.

Y. Chen, S. Bhojanapalli, S. Sanghavi, and R. Ward. Coherent matrix completion. In Proceedings of the 31st International Conference on Machine Learning, pages 674–682. JMLR, 2014.

Y. Chen, Y. Chi, J. Fan, C. Ma, and Y. Yan. Noisy matrix completion: Understanding statistical guarantees for convex relaxation via nonconvex optimization. arXiv preprint arXiv:1902.07698, 2019b.

Figure 11: Coherent and Nonuniform Sampling (NUS) simulations with SNR = 10. The optimal nonconvex penalties are both obtained at = 5 under the two scenarios respectively. The integers from 1 to 100 on the x-axis index the grid of 100 values of (from largest to smallest) as described in Section 4.1.

Y. Chi, Y. M. Lu, and Y. Chen. Nonconvex optimization meets low-rank matrix factorization: An overview. IEEE Transactions on Signal Processing, 67(20):5239–5269, 2019.

A. L. Chistov and D. Y. Grigor’ev. Complexity of quantifier elimination in the theory of algebraically closed fields. In Mathematical Foundations of Computer Science 1984, pages 17–31. Springer, 1984.

I. Daubechies, R. DeVore, M. Fornasier, and C. S. G¨unt¨urk. Iteratively reweighted least squares minimization for sparse recovery. Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 63(1):1–38, 2010.

A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society, Series B, 39:1–38, 1977.

B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least angle regression (with discussion). Annals of Statistics, 32(2):407–499, 2004. ISSN 0090-5364.

(a) ROM, 90% missing, SNR = 1, true rank = 10

(d) ROM, 95% missing, SNR = 5, true rank = 5

Figure 12: The y-axis denotes the number of iterations NC-Impute takes to stabilize the rank. The integers on the x-axis index some values on a grid of (from largest to smallest) as described in Section 4.1. The six plots represent the six scenarios considered in Section 4.1: (a)-(d) correspond to the four scenarios of Example-A; (e) covers Example-B; (f) is for Example-C. Each procedure is repeated 10 times.

J. Fan and R. Li. Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association, 96:1348–1360, 2001.

M. Fazel. Matrix Rank Minimization with Applications. PhD thesis, Stanford University, 2002.

L. Feng and C.-H. Zhang. Sorted concave penalized regression. arXiv preprint arXiv:1712.09941, 2017.

M. Fornasier, H. Rauhut, and R. Ward. Low-rank matrix recovery via iteratively reweighted least squares minimization. SIAM Journal on Optimization, 21(4):1614–1640, 2011.

I. E. Frank and J. H. Friedman. A statistical view of some chemometrics regression tools. Technometrics, 35:109–135, 1993.

R. M. Freund, P. Grigas, and R. Mazumder. An Extended Frank-Wolfe Method with “In-Face” Directions, and its Application to Low-Rank Matrix Completion. ArXiv e-prints, November 2015.

R. Ge, J. D. Lee, and T. Ma. Matrix completion has no spurious local minimum. In Advances in Neural Information Processing Systems, pages 2973–2981, 2016.

R. Ge, C. Jin, and Y. Zheng. No spurious local minima in nonconvex low rank problems: A unified geometric analysis. In Proceedings of the 34th International Conference on Machine LearningVolume 70, pages 1233–1242. JMLR. org, 2017.

G. Golub and C. Van Loan. Matrix Computations. Johns Hopkins University Press, Baltimore., 1983.

D. Gross. Recovering low-rank matrices from few coefficients in any basis. IEEE Transactions on Information Theory, 57(3):1548–1566, 2011.

S. Gu, Q. Xie, D. Meng, W. Zuo, X. Feng, and L. Zhang. Weighted nuclear norm minimization and its applications to low level vision. International journal of computer vision, 121(2):183–208, 2017.

M. Hardt. Understanding alternating minimization for matrix completion. In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pages 651–660. IEEE, 2014.

M. Hardt and M. Wootters. Fast matrix completion without the condition number. In Conference on learning theory, pages 638–678, 2014.

T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning: Prediction, Inference and Data Mining (Second Edition). Springer Verlag, New York, 2009.

T. Hastie, R. Mazumder, J. D. Lee, and R. Zadeh. Matrix completion and low-rank svd via fast alternating least squares. Journal of Machine Learning Research, to appear, 2016.

H. Hazimeh and R. Mazumder. Fast best subset selection: Coordinate descent and local combinatorial optimization algorithms. Operations Research, to appear, 2019.

R. A. Horn and C. R. Johnson. Matrix analysis. Cambridge university press, 2012.

M. Jaggi and M. Sulovsk´y. A simple algorithm for nuclear norm regularized problems. In Proceedings of the 27th International Conference on Machine Learning (ICML-10), pages 471–478, 2010.

P. Jain, R. Meka, and I. S. Dhillon. Guaranteed rank minimization via singular value projection. In Advances in Neural Information Processing Systems, pages 937–945, 2010.

P. Jain, P. Netrapalli, and S. Sanghavi. Low-rank matrix completion using alternating minimization. In Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pages 665–674. ACM, 2013.

R. H. Keshavan, A. Montanari, and S. Oh. Matrix completion from noisy entries. Journal of Machine Learning Research, 11:2057–2078, 2010.

O. Klopp. Noisy low-rank matrix completion with general sampling distribution. Bernoulli, 20(1): 282–303, 2014.

V. Koltchinskii, K. Lounici, and A. B. Tsybakov. Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion. The Annals of Statistics, 39(5):2302–2329, 2011.

R. Larsen. Propack-software for large and sparse svd calculations, 2004. Available at

G. Lecu´e and S. Mendelson. Regularization and the small-ball method i: sparse recovery. The Annals of Statistics, 46(2):611–641, 2018.

A. S. Lewis. The convex analysis of unitarily invariant matrix functions. Journal of Convex Analysis, 2:173–183, 1995.

P.-L. Loh and M. J. Wainwright. Regularized m-estimators with nonconvexity: Statistical and algo- rithmic theory for local optima. Journal of Machine Learning Research, 16:559–616, 2015.

J. Lv and Y. Fan. A unified approach to model selection and sparse recovery using regularized least squares. The Annals of Statistics, 37:3498–3528, 2009.

C. Ma, K. Wang, Y. Chi, and Y. Chen. Implicit regularization in nonconvex statistical estimation: Gradient descent converges linearly for phase retrieval, matrix completion and blind deconvolution. arXiv preprint arXiv:1711.10467, 2017.

R. Mazumder and P. Radchenko. The discrete dantzig selector: Estimating sparse linear models via mixed integer linear optimization. arXiv preprint arXiv:1508.01922, 2015.

R. Mazumder and H. Weng. Computing the degrees of freedom of rank-regularized estimators and cousins. Electronic Journal of Statistics, to appear, 2020.

R. Mazumder, T. Hastie, and R. Tibshirani. Spectral regularization algorithms for learning large incomplete matrices. Journal of Machine Learning Research, 11:2287–2322, 2010.

R. Mazumder, J. H. Friedman, and T. Hastie. Sparsenet: coordinate descent with nonconvex penalties. Journal of the American Statistical Association, 106:1125–1138, 2011.

R. Mazumder, P. Radchenko, and A. Dedieu. Subset selection with shrinkage: Sparse linear modeling when the snr is low. arXiv preprint arXiv:1708.03288, 2017.

K. Mohan and M. Fazel. Reweighted nuclear norm minimization with application to system identifi- cation. In Proceedings of the 2010 American Control Conference, pages 2953–2959. IEEE, 2010.

K. Mohan and M. Fazel. Iterative reweighted algorithms for matrix rank minimization. Journal of Machine Learning Research, 13(Nov):3441–3473, 2012.

S. N. Negahban and M. J. Wainwright. Estimation of (near) low-rank matrices with noise and high- dimensional scaling. The Annals of Statistics, 39:1069–1097, 2011.

S. N. Negahban and M. J. Wainwright. Restricted strong convexity and weighted matrix completion: optimal bounds with noise. Journal of Machine Learning Research, 13:1665–1697, 2012.

M. Nikolova. Local strong homogeneity of a regularized estimator. SIAM Journal on Applied Mathematics, 61:633–658, 2000.

B. Recht. A simpler approach to matrix completion. Journal of Machine Learning Research, 12: 3413–3430, 2011.

B. Recht, M. Fazel, and P. A. Parrilo. Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM Review, 52:471–501, 2010.

R. T. Rockafellar. Convex Analysis. Princeton University Press, Princeton, New Jersey, 1970.

A. Rohde and A. B. Tsybakov. Estimation of high-dimensional low-rank matrices. The Annals of Statistics, 39:887–930, 2011.

A. Shapiro, Y. Xie, and R. Zhang. Matrix completion with deterministic pattern: A geometric perspective. IEEE Transactions on Signal Processing, 67(4):1088–1103, 2018.

A. SIGKDD and Netflix. Soft modelling by latent variables: the nonlinear iterative partial least squares (NIPALS) approach. In Proceedings of KDD Cup and Workshop, 2007.

C. M. Stein. Estimation of the mean of a multivariate normal distribution. The Annals of Statistics, pages 1135–1151, 1981.

G. W. Stewart and J.-G. Sun. Matrix Perturbation Theory. Computer science and scientific computing. Academic Press, 1990.

R. Sun and Z.-Q. Luo. Guaranteed matrix completion via non-convex factorization. IEEE Transactions on Information Theory, 62(11):6535–6579, 2016.

R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B, 58:267–288, 1996.

O. Troyanskaya, M. Cantor, G. Sherlock, P. Brown, T. Hastie, R. Tibshirani, D. Botstein, and R. B. Altman. Missing value estimation methods for dna microarrays. Bioinformatics, 17(6):520–525, 2001.

S. Wang, H. Weng, and A. Maleki. Which bridge estimator is optimal for variable selection? Annals of Statistics, to appear, 2019.

H. Weng, A. Maleki, L. Zheng, et al. Overcoming the limitations of phase transition by higher order analysis of regularization techniques. The Annals of Statistics, 46(6A):3099–3129, 2018.

C.-H. Zhang. Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics, 38:894–942, 2010.

C.-H. Zhang and T. Zhang. A general theory of concave regularization for high-dimensional sparse estimation problems. Statistical Science, 27(4):576–593, 2012.

L. Zheng, A. Maleki, H. Weng, X. Wang, and T. Long. Does -minimization outperform minimization? IEEE Transactions on Information Theory, 63(11):6896–6935, 2017.

Q. Zheng and J. Lafferty. Convergence analysis for rectangular matrix completion using burer-monteiro factorization and gradient descent. arXiv preprint arXiv:1605.07051, 2016.

H. Zou. The adaptive lasso and its oracle properties. Journal of the American statistical association, 101(476):1418–1429, 2006.

H. Zou and R. Li. One-step sparse estimates in nonconcave penalized likelihood models. The Annals of Statistics, 36(4):1509–1533, 2008.