Student, IEEE Mohammadreza Malek-Mohammadi

Senior Member, IEEE Massoud Babaie-Zadeh

Arash Amini
email: aamini@sharif.edu

Fellow, IEEE Christian Jutten
email: christian.jutten@inpg.fr

Recovery of Low Rank Matrices Under Affine Constraints via a Smoothed Rank Function

Keywords: Affine Rank Minimization (ARM), Compressive Sensing, Matrix Completion (MC), Nuclear Norm Minimization (NNM), Rank Approximation, Spherical Section Property (SSP). I

In this paper, the problem of matrix rank minimization under affine constraints is addressed. The state-of-theart algorithms can recover matrices with a rank much less than what is sufficient for the uniqueness of the solution of this optimization problem. We propose an algorithm based on a smooth approximation of the rank function, which practically improves recovery limits on the rank of the solution. This approximation leads to a non-convex program; thus, to avoid getting trapped in local solutions, we use the following scheme. Initially, a rough approximation of the rank function subject to the affine constraints is optimized. As the algorithm proceeds, until reaching the desired accuracy, finer approximations of the rank are successively optimized while the solver is initialized with the solution of the previous approximation.

On the theoretical side, benefiting from the spherical section property, we will show that the sequence of the solutions of the approximating programs converges to the minimum rank solution. On the experimental side, it will be shown that the proposed algorithm, termed SRF standing for Smoothed Rank Function, can recover matrices which are unique solutions of the rank minimization problem and yet not recoverable by nuclear norm minimization. Furthermore, it will be demonstrated that, in completing partially observed matrices, the accuracy of SRF is considerably and consistently better than some famous algorithms when the number of revealed entries is close to the minimum number of parameters that uniquely represent a lowrank matrix.

certain random laws and the rank of the matrix is sufficiently small, one can uniquely recover the matrix with overwhelming probability [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]- [START_REF] Recht | Null space conditions and thresholds for rank minimization[END_REF].

One of the celebrated applications of affine rank minimization (or matrix completion) is Collaborative Filtering [START_REF] Candés | Exact matrix completion via convex optimization[END_REF]. This technique is applied when a system tries to recommend goods to customers/users based on the available feedbacks of all the customers. In fact, the system learns the user preferences through the feedbacks and identifies similarities between them. As the number of factors affecting the user interests is much less than the total number of customers and products, the matrix whose (i, j)-th entry represents the rating of the ith user for the j-th product is expected to be low-rank. This could be efficiently used by the matrix completion techniques to predict the users' ratings for unrated items.

Applications of affine rank minimization in control theory include System Identification [START_REF] Candés | Matrix completion with noise[END_REF] and low-order realization of linear systems [START_REF] Fazel | A rank minimization heuristic with application to minimum order system approximation[END_REF]. In the former, the goal is to find an LTI system with minimum order that fits the available joint input-output observations of a multiple-input multiple-output system [START_REF] Zh | Interior-point method for nuclear norm approximation with application to system identification[END_REF].

In wireless sensor networks, due to limited energy resources and transmitting power, the sensors are able to communicate only with their neighboring sensors. These communications (e.g., received powers) determine the pairwise distances between sensors, which partially reveals the matrix of all pairwise distances. To localize the sensors in the network, one needs to estimate their distances from predefined anchor points which in turn requires completion of the distance matrix through the multi-dimensional scaling technique [START_REF] Cox | Multidimensional Scaling[END_REF]. Interestingly, the rank of the pairwise distance matrix is small compared to its dimension [START_REF] Candés | Matrix completion with noise[END_REF].

Other areas to which affine rank minimization is applied include Machine Learning [START_REF] Amit | Uncovering shared structures in multiclass classification[END_REF], Quantum State Tomography [START_REF] Gross | Quantum state tomography via compressed sensing[END_REF], Spectrum Sensing [START_REF] Meng | Collaborative spectrum sensing from sparse observations in cognitive radio networks[END_REF], and Spatial Covariance matrix completion [START_REF] Candés | Matrix completion with noise[END_REF], [START_REF] Ito | Crystalmusic: Accurate localization of multiple sources in diffuse noise environments using crystal-shaped microphone arrays[END_REF]. The spatial covariance matrix is essential in estimating the directions of arrival of sources impinging on an array of antennas using for example MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] or ESPRIT algorithms [START_REF] Roy | Esprit-estimation of signal parameters via rotational invariance techniques[END_REF].

The main difficulty of the affine rank minimization problem is due to the fact that the rank function is discontinuous and non-differentiable. Indeed, the optimization problem is NPhard, and all available optimizers have doubly exponential complexity [START_REF] Chistov | Complexity of quantifier elimination in the theory of algebraically closed fields[END_REF]. In [START_REF] Fazel | Matrix Rank Minimization with Applications[END_REF], Fazel proposed to replace the rank of the matrix with its nuclear norm, which is defined as the sum of all singular values (SV). This modification is known to be the tightest convex relaxation of the rank minimization problem [START_REF] Candés | Matrix completion with noise[END_REF] and can be implemented using a Semi Definite Program (SDP) [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]. Using similar techniques as in compressed sensing, it is recently shown that under mild conditions and with overwhelming probability, the nuclear norm minimization (NNM) technique achieves the same solution as the original rank minimization approach [START_REF] Mohan | A simplified approach to recovery conditions for low rank matrices[END_REF]- [START_REF] Kong | Sufficient conditions for low-rank matrix recovery, translated from sparse signal recovery[END_REF].

Other approaches toward rank minimization consist of either alternative solvers instead of SDP in NNM or approximating the rank function using other forms rather than the nuclear norm. The FPCA method belongs to the first category and uses fixed point and Bergman iterative algorithm to solve NNM [START_REF] Ma | Fixed point and bregman iterative methods for matrix rank minimization[END_REF]. Among the examples of the second category, one can name LMaFit [START_REF] Wen | Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm[END_REF], BiG-AMP [START_REF] Parker | Bilinear generalized approximate message passing[END_REF], and OptSpace [START_REF] Keshavan | Matrix completion from a few entries[END_REF]. It is also possible to generalize the greedy methods of compressive sensing to the rank minimization problem; for instance, ADMiRA [START_REF] Lee | Admira: Atomic decomposition for minimum rank approximation[END_REF] generalizes the CoSaMP [START_REF] Needell | Cosamp: Iterative signal recovery from incomplete and inaccurate samples[END_REF].

In this work, we introduce an iterative method that is based on approximating the rank function. However, in contrast to previous methods, the approximation is continuous and differentiable, is made finer at each iteration, and, asymptotically, will coincide with the rank function. Our method is inspired by the work of Mohimani et al [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF] which uses smoothed ℓ 0norm1 to obtain sparse solutions of underdetermined system of linear equations. Nevertheless, the way SRF is extended from [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF], and, particularly, the performance guarantees that are provided are among the contribution of our paper. Furthermore, in generalizing the method of [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF] to the ARM problem, we need to derive the gradient of the rank approximating functions in a closed form which is another novelty of the current work.

A few preliminary results of this work have been presented in the conference paper [START_REF] Ghasemi | SRF: Matrix completion based on smoothed rank function[END_REF]. While [START_REF] Ghasemi | SRF: Matrix completion based on smoothed rank function[END_REF] was only devoted to the matrix completion problem, the current paper focuses on the more general problem of affine rank minimization. Furthermore, here, we present mathematical and experimental convergence analysis, and consider more comprehensive numerical evaluation scenarios.

The reminder of this paper is organized as follows. In Section II, the ARM problem is formulated, and in Section III, the SRF algorithm is introduced. Section IV is devoted to analyze the convergence properties of the SRF algorithm. In Section V, some experimental results of our algorithm are provided, and it will be compared empirically against some well known algorithms. Finally, Section VI concludes the paper.

II. PROBLEM FORMULATION

The affine rank minimization problem generally is formulated as min

X rank(X) subject to A(X) = b, (1)
where X ∈ R n1×n2 is the decision variable, A : R n1×n2 → R m is a known linear operator, and b ∈ R m is the observed measurement vector. The affine constraints A(X) = b can be converted to

A vec(X) = b, (2)
where A ∈ R m×n1n2 denotes the matrix representation of the linear operator A and vec(X) denotes the vector in R n1n2 with the columns of X stacked on top of one another. The special case of matrix completion corresponds to the setting

min X rank(X) subject to [X] ij = [M] ij ∀(i, j) ∈ Ω, (3)
where X is as in [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF], M ∈ R n1×n2 is the matrix whose entries are partially observed, Ω ⊂ {1, 2, ..., n 1 } × {1, 2, ..., n 2 } is the set of the indexes of the observed entries of M, and [X] ij is the (i, j)-th entry of X. Indeed, the constraints [X] ij = [M] ij ,∀(i, j) ∈ Ω is an affine mapping which keeps some of the entries and discards others.

In the nuclear norm minimization, the rank function is replaced with the nuclear norm of the decision variable, leading to

min X X * subject to A(X) = b, (4)
where X * r i=1 σ i (X) is the nuclear norm, in which r is the rank of the matrix X, and σ i (X) is the i-th largest singular value of the matrix X. There is a strong parallelism between this rank minimization and ℓ 0 -norm minimization in compressive sensing [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]. In particular, minimizing the rank is equivalent to minimizing the number of non-zero singular values. Hence, (1) can be reformulated as

min X σ(X) 0 subject to A(X) = b, (5)
where σ(X) = (σ 1 (X), ..., σ n (X)) T is the vector of all singular values, • 0 denotes the ℓ 0 -norm, and n = min(n 1 , n 2 ). 2Likewise, the nuclear norm is the ℓ 1 -norm of the singular value vector where ℓ 1 -norm of a vector, denoted by • 1 , is the sum of the absolute values of its elements. This suggests the alternative form of

min X σ(X) 1 subject to A(X) = b (6)
for [START_REF] Candés | Matrix completion with noise[END_REF]. Based on this strong parallel, many results in compressive sensing theory (see for example [START_REF] Candès | Decoding by linear programming[END_REF]- [START_REF] Zhang | Theory of compressive sensing via ℓ 1 minimization: A nonrip analysis and extensions[END_REF]) have been adopted in the rank minimization problem [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF], [START_REF] Mohan | A simplified approach to recovery conditions for low rank matrices[END_REF], [START_REF] Oymak | New null space results and recovery thresholds for matrix rank minimization[END_REF], [START_REF] Oymak | Improved thresholds for rank minimizations[END_REF].

III. THE PROPOSED ALGORITHM

A. The main idea

Our approach to solve the ARM problem is to approximate the rank with a continuous and differentiable function, and then to use a gradient descent algorithm to minimize it. The approximation is such that the error can be made arbitrarily small. In contrast, note that the nuclear norm is not differentiable [START_REF] Lewis | Nonsmooth analysis of singular values. part i: Theory[END_REF] and its approximation error depends on the singular values of the matrix and cannot be controlled.

Instead of using a fixed approximation, we use a family G δ : R n1×n2 → R + of approximations, where the index δ is a measure of approximation error and reflects the accuracy. The smaller δ, the closer behavior of G δ to the rank. For instance, G 0 stands for the errorless approximation; i.e., G 0 coincides with the rank function. We constrain the family to be continuous with respect to δ. This helps in achieving the rank minimizer (G 0 ) by gradually decreasing δ. Besides, to facilitate finding the minimizers of the relaxed problem, we require the G δ 's for δ > 0 to be differentiable with respect to the input matrix.

In order to introduce suitable G δ families, we specify certain families of one-dimensional functions that approximate Kronecker delta function.

Assumption 1: Let f : R → [0, 1] and define f δ (x) = f (x/δ) for all δ > 0. The class {f δ } is said to satisfy the Assumption 1, if (a) f is real, symmetric unimodal, and analytic, (b)

f (x) = 1 ⇔ x = 0, (c) f ′′ (0) < 0, and (d) lim |x|→∞ f (x) = 0.
It follows from Assumption 1 that x = 0 is the unique mode of all f δ 's. This implies that f ′ δ (0) = 0 for δ = 0. In addition, {f δ } converge pointwise to Kronecker delta function as δ → 0, i.e.,

lim δ→0 f δ (x) = 0 if x = 0, 1 if x = 0. (7)
The class of Gaussian functions, which is of special interest in this paper, is defined as

f δ (x) = exp(- x 2 2δ 2 ). (8)
It is not difficult to verify the constraints of Assumption 1 for this class. Other examples include

f δ (x) = 1 -tanh( x 2 2δ 2 ) and f δ (x) = δ 2
x 2 +δ 2 . To extend the domain of {f δ } to matrices, let define

F δ (X) = h δ σ(X) = n i=1 f δ σ i (X) , (9)
where n = min(n 1 , n 2 ) and h δ : R n → R is defined as

h δ (x) = n i=1 f δ (x i ).
Since f δ is an approximate Kronecker delta function, F δ (X) yields an estimate of the number of zero singular values of X. Consequently, it can be concluded that rank(X) ≈ n -F δ (X), and the ARM problem can be relaxed to min

X G δ (X) = n -F δ (X) subject to A(X) = b, (10) or equivalently max X F δ (X) subject to A(X) = b. ( 11
)
The advantage of maximizing F δ compared to minimizing the rank is that F δ is smooth and we can apply gradient methods. However, for small values of δ where G δ is a relatively good approximate of the rank function, F δ has many local maxima, which are likely to trap gradient methods.

To avoid local maxima 3 , we initially apply a large δ. Indeed, we will show in Theorem 2 that under, Assumption 1, F δ becomes concave as δ → ∞ and (11) will have a unique solution. Then we gradually decrease δ to improve the accuracy of approximation. For each new value of δ, we initialize the maximization of F δ with the result of [START_REF] Ito | Crystalmusic: Accurate localization of multiple sources in diffuse noise environments using crystal-shaped microphone arrays[END_REF] for the previous value of δ. From the continuity of {f δ } with respect to δ, it is expected that the solutions of [START_REF] Ito | Crystalmusic: Accurate localization of multiple sources in diffuse noise environments using crystal-shaped microphone arrays[END_REF] for δ i and δ i+1 are close, when δ i and δ i+1 are close. In this fashion, the chance of finding a local maximum instead of a global one is decreased. This approach for optimizing non-convex functions is known as Graduated Non-Convexity (GNC) [START_REF] Blake | Visual Reconstruction[END_REF], and was used in [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF] to minimize functions approximating the ℓ 0 -norm.

B. Gradient Projection

For each δ in the decreasing sequence, to maximize F δ with equality constraints, we use the Gradient Projection (GP) technique [START_REF] Bertsekas | Nonlinear Programming[END_REF]. In GP, the search path at each iteration is obtained by projecting back the ascent (or descent) direction onto the feasible set [START_REF] Bertsekas | Nonlinear Programming[END_REF]. In other words, at each iteration, one has X ← P X + µ j ∇F δ (X) , where P denotes the orthogonal projection onto the affine set defined by linear constraints A(X) = b, and µ j is the step-size of the j-th iteration. As the feasible set is affine, several methods can be exploited to implement the projection P. For example, one can store the QR factorization of the matrix implementation of A for fast implementation of the back projection, or, alternatively, a least-squares problem can be solved at each step [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]. The closed form solution of the least-squares problem can be found in Appendix A.

To complete the GP step, we should derive the gradient of the approximating functions with respect to the matrix X. Surprisingly, although σ i (X), i = 1, ..., n and X * are not differentiable functions of X [START_REF] Lewis | Nonsmooth analysis of singular values. part i: Theory[END_REF], the following theorem shows that one can find functions F δ = h δ • σ(X) which are differentiable under the absolutely symmetricity of the h δ . Before stating the theorem, recall that a function f : R q → [-∞, +∞] is called absolutely symmetric [START_REF] Lewis | The convex analysis of unitarily invariant matrix norms[END_REF] if f (x) is invariant under arbitrary permutations and sign changes of the components of x.

Theorem 1: Suppose that

F δ : R n1×n2 → R is rep- resented as F δ (X) = h δ σ(X) = h δ • σ(X), where X ∈ R n1×n2 with the Singular Value Decomposition (SVD) X = Udiag(σ 1 , ..., σ n )V T , σ(X) : R n1×n2 → R n has the SVs of the matrix X, n = min(n 1 , n 2 )
, and h δ : R n → R is absolutely symmetric and differentiable. Then the gradient of

F δ (X) at X is ∂F δ (X) ∂X = Udiag(θ)V T , (12)
where θ = ∂h δ (y) ∂y | y=σ(X) denotes the gradient of h δ at σ(X). Informal Proof : In [START_REF] Lewis | The convex analysis of unitarily invariant matrix norms[END_REF]Cor. 2.5], it is shown that if a function h δ is absolutely symmetric and the matrix X has σ(X) in the domain of h δ , then the subdifferential4 of F δ is given by

∂ h δ • σ(X) = {Udiag(θ)V T |θ ∈ ∂h δ σ(X) }. (13)
Since h δ is differentiable, ∂h δ σ(X) is a singleton and consequently ∂ h δ • σ(X) becomes a singleton. When the subdifferential of a non-convex function becomes singleton, the function is intuitively expected to be differentiable with the subgradient as its gradient. 5 Nevertheless, to the best of our knowledge, there is no formal proof. Provided that this intuition is true, then ∂ h δ • σ(X) will be converted to ∇ h δ • σ(X) and equation ( 12) is obtained.

Formal Proof: Equation ( 12) can be obtained directly from the "if part" of [START_REF] Lewis | The convex analysis of unitarily invariant matrix norms[END_REF]Thm. 3.1], which does not require convexity of h δ as stated in its proof.

Corollary 1: For the Gaussian function family given in ( 8), the gradient of

F δ (X) at X is ∂F δ (X) ∂X = Udiag(- σ 1 δ 2 e -σ 2 1 /2δ 2 , ..., - σ n δ 2 e -σ 2 n /2δ 2 )V T . ( 14
)
Proof: f δ is an even function for the Gaussian family; therefore, h δ becomes an absolutely symmetric function. As a result, Theorem 1 proves [START_REF] Chistov | Complexity of quantifier elimination in the theory of algebraically closed fields[END_REF].

C. Initialization

Naturally, we initialize the GNC procedure by the solution of ( 11) corresponding to δ → ∞. This solution can be found from the following theorem.

Theorem 2: Consider a class of one variable functions {f δ } satisfying the Assumption 1. For the rank approximation problem [START_REF] Ito | Crystalmusic: Accurate localization of multiple sources in diffuse noise environments using crystal-shaped microphone arrays[END_REF]

, let X = argmin{ X F | A(X) = b}, then lim δ→∞ argmax{F δ (X) | A(X) = b} = X, (15)
where • F denotes the matrix Frobenius norm.

There is a simple interpretation of the solution of (11) for the Gaussian family when δ approaches ∞. As e -x ≈ 1 -x for small values of x,

F δ (X) = n i=1 e -σ 2 i (X)/2δ 2 ≈ n - n i=1 σ 2 i (X)/δ 2 for δ ≫ σ i (X). Consequently, argmax{F δ (X) | A(X) = b} ≈ argmin{ n i=1 σ 2 i (X)|A(X) = b} = argmin{ X F |A(X) = b}.
The proof is left to Appendix B.

The following Corollary is an immediate result of the above theorem.

Corollary 2: For the matrix completion problem, the initial solution of the SRF algorithm is X with the following definition:

[ X] ij = [M] ij (i, j) ∈ Ω, 0 (i, j) / ∈ Ω. , (16)
where M and Ω are as defined in (3).

• Initialization: 1) Let X0 = argmin{ X F | A(X) = b} as the initial solution. 2) Choose a suitable decreasing sequence of δ, {δ1, δ2, . . .}; e.g., δj = cδj-1, j ≥ 2. 3) Choose ǫ as the stopping threshold. 4) Choose suitable L (Number of internal loop iteration) and µ, and initialize j with 1.

• While d > ǫ 1) Let δ = δj.

2) Internal maximization loop:

-Initialization: X = Xj-1.

-For ℓ = 1 . . . L, a) Compute the SVD of

X = Udiag(σ1, . . . , σn)V T . b) Let D = Udiag(-σ1e -σ 2 1 /2δ 2 , . . . , -σne -σ 2 n /2δ 2 )V T . c) X ← X + µD.
d) Project X back onto the feasible set:

X ← P(X). 3) Set Xj = X. 4) d = Xj -Xj-1 F / √ n1n2. 5) j ← j + 1.
• Final answer is Xj.

D. The Final Algorithm

The final algorithm is obtained by applying the main idea, initial solution, and gradient projection to the Gaussian function given in [START_REF] Amit | Uncovering shared structures in multiclass classification[END_REF]. Fig. 1 depicts the algorithm. In the sequel, we briefly review some remarks about the parameters used in the implementation of the algorithm. Most of these remarks correspond to similar remarks for the SL0 algorithm [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF] and are presented here for the sake of completeness.

Remark 1. It is not necessary to wait for the convergence of the internal steepest ascent loop because as explained in Section III-A for each value of δ, it is just needed to get close to the global maximizer of F δ to avoid local maxima. Therefore, the internal loop is only repeated for a fixed number of times (L).

Remark 2. After initiating the algorithm with the minimum Frobenius norm solution, the first value of δ may be set to about two to four times of the largest SV of X 0 (the initial guess). If we take δ > 4 max i σ i ( X 0 ) , then exp -σ 2 i ( X 0 )/2δ 2 > 0.97 ≈ 1 for 1 ≤ i ≤ n. Thus, this δ value acts virtually like ∞ for all SVs of X 0 . In addition, the decreasing sequence can be adjusted to δ j = cδ j-1 , j ≥ 2, where c generally is chosen between 0.5 and 1.

Remark 3. This remark is devoted to the selection of µ j , step-size parameter. Typically, in a gradient ascent algorithm, µ j should be chosen small enough to follow the ascent direction. Furthermore, reducing δ results in more fluctuating behaviour of the rank approximating function. Therefore, to avoid large steps which cause jumps over the maximizer, one should choose smaller values of step-size for smaller values of δ. Following the same reasoning as in [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF]Remark 2], a good choice is to decrease µ j proportional to δ 2 ; that is, µ j = µδ 2 , where µ is a constant. sBy letting µ j = µδ 2 , the gradient step can be reduced to

X ← X -µUdiag(σ 1 e -σ 2 1 /2δ 2 , . . . , σ n e -σ 2 n /2δ 2 )V T .
Remark 4. The distance between the solutions at the two consecutive iterations is the criterion to stop the algorithm. That is, if d X j -X j-1 F / √ n 1 n 2 is smaller than some tolerance (ǫ), the iterations are ended and X j becomes the final solution.

IV. CONVERGENCE ANALYSIS

Noting that the original problem is NP-hard and we are dealing with maximizing non-concave functions, a complete and thorough convergence analysis would be beyond the scope of this paper. We believe that similar to [START_REF] Mohimani | Sparse recovery using smoothed ℓ 0 (SL0): Convergence analysis[END_REF] which examines the global convergence properties of the SL0 algorithm [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF], it would be possible to analyze the convergence of the SRF algorithm to the global solution. However, in this paper, we only study a simplified convergence analysis, and the complete analysis is left for a future work.

For the simplified analysis, in the sequel, it is assumed that the internal loop has been converged to the global maximum, and we prove that this global solution converges to the minimum rank solution as δ goes to zero. This analysis helps us to characterize the conditions under which

lim δ→0 argmax{F δ (X) | A(X) = b} (17)
is equivalent to

argmax{lim δ→0 F δ (X) | A(X) = b}. ( 18
)
The equivalence of ( 17) and ( 18) is of particular importance since it shows that the idea of SRF corresponding to optimization of ( 17) is indeed the case and leads to finding the solution of program [START_REF] Kong | Sufficient conditions for low-rank matrix recovery, translated from sparse signal recovery[END_REF] which is identical to the original affine rank minimization problem defined in [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF].

The following results and proofs are not direct extension of the convergence results of [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF] and are more tricky to obtain, though our exposition follows the same line of presentation.

We start the convergence analysis by the definition of Spherical Section Property (SSP), used in the analysis of uniqueness of the rank and nuclear norm minimization [START_REF] Mohan | A simplified approach to recovery conditions for low rank matrices[END_REF], and a lemma which makes this abstract definition clearer.

Definition 1: Spherical Section Property [START_REF] Mohan | A simplified approach to recovery conditions for low rank matrices[END_REF], [START_REF] Dvijotham | A nullspace analysis of the nuclear norm heuristic for rank minimization[END_REF]. The spherical section constant of a linear operator A : R n1×n2 → R m is defined as

∆(A) = min Z∈null(A)\{0} Z 2 * Z 2 F . (19)
Further, A is said to have the ∆-spherical section property if

∆(A) ≥ ∆.
Definition 1 extends a similar concept in the compressive sensing framework where it is shown that many randomly generated sensing matrices possesses SSP with high probability [START_REF] Zhang | Theory of compressive sensing via ℓ 1 minimization: A nonrip analysis and extensions[END_REF]. Although extending a similar theoretical result to the matrix case is a topic of interest, [START_REF] Dvijotham | A nullspace analysis of the nuclear norm heuristic for rank minimization[END_REF] proves that if all entries of the matrix representation of A are identically and independently distributed (i.i.d.) from a zero-mean, unit-variance Gaussian distribution, then, under some mild conditions, A possesses the ∆-spherical section property with overwhelming probability.

Lemma 1: Assume A has the ∆-spherical section property. Then for any X ∈ null(A) \ {0}, we have rank(X) ≥ ∆.

Proof: Since X belongs to null(A), one can write

X * X F ≥ √ ∆ ⇒ X * ≥ √ ∆ X F .
It is also known that rank(X) X F ≥ X * , see for example [START_REF] Horn | Matrix analysis[END_REF]. Putting them together, we have

X * ≥ √ ∆ X * rank(X) ⇒ rank(X) ≥ ∆
or rank(X) ≥ ⌈∆⌉, where ⌈∆⌉ denotes the smallest integer greater than or equal to ∆.

The above lemma shows that if ∆ is large, the null space of A does not include low rank matrices. Such subspaces are also known as almost Euclidean subspaces [START_REF] Zhang | Theory of compressive sensing via ℓ 1 minimization: A nonrip analysis and extensions[END_REF], in which the ratio of ℓ 1 -norm to ℓ 2 -norm of elements cannot be small. Theorem 3 ([39]): Suppose A has the ∆-spherical property, and

X 0 ∈ R n1×n2 satisfies A(X 0 ) = b. If rank(X 0 ) < ∆ 2
, then X 0 is the unique solution of problem [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF].

Lemma 2: Assume A : R n1×n2 → R m has the ∆-spherical section property, and set n = min(n 1 , n 2 ). Let X be any element in null(A) and σ 1 , ..., σ n represent its singular values. Then for any subset I of {1, ..., n} such that |I| + ∆ > n,

i∈I σ i ( n i=1 σ 2 i ) 0.5 ≥ √ ∆ -n -|I|, (20)
where | • | denotes the cardinality of a set. Proof:

If I = {1, ..., n}, then it is clear that n i=1 σi ( n i=1 σ 2 i ) 0.5 ≥ √ ∆, since the ∆-spherical section property holds. Otherwise, if |I| < n, the ∆-spherical section property implies that √ ∆ ≤ X * X F = n i=1 σ i ( n i=1 σ 2 i ) 0.5 .
For the sake of simplicity, let us define

α i = σ i ( n i=1 σ 2 i ) 0.5 . This shows that 1 = n i=1 α 2 i ≥ i / ∈I α 2 i ≥ ( i / ∈I α i ) 2 n -|I| ,
where we used the inequality ∀z ∈ R p , z 2 1 ≤ p z 2 2 [START_REF] Horn | Matrix analysis[END_REF]. Hence, it can be concluded that

i / ∈I α i ≤ n -|I|.
On the other hand, it is known that

√ ∆ ≤ i∈I α i + i / ∈I α i ≤ i∈I α i + n -|I|, which confirms that i∈I σ i ( n i=1 σ 2 i ) 0.5 = i∈I α i ≥ √ ∆ -n -|I|.
Corollary 3: If A : R n1×n2 → R m has the ∆-spherical section property, n = min(n 1 , n 2 ), and X ∈ null(A) has at most ⌈∆ -1⌉ singular values greater than α, then

X F ≤ nα √ ∆ -⌈∆ -1⌉ .
Proof: At least n -⌈∆ -1⌉ singular values of X are less than or equal to α. If I denotes the indices of singular values not greater than α, then by using Lemma 2, we will have

i∈I σ i ( n i=1 σ 2 i ) 0.5 ≥ √ ∆ -n -n + ⌈∆ -1⌉ ⇒ X F ( √ ∆ -⌈∆ -1⌉) ≤ i∈I σ i ≤ nα,
which proves that

X F ≤ nα √ ∆ -⌈∆ -1⌉ .
Lemma 3: Assume A : R n1×n2 → R m has the ∆-spherical section property, f δ (•) is a member of the class that satisfies Assumption 1, and define F δ as in [START_REF] Gross | Quantum state tomography via compressed sensing[END_REF] and n = min(n 1 , n 2 ). Let X = {X|A(X) = b} contain a solution X 0 with rank(X 0 ) = r 0 < ∆ 2 . Then for any X ∈ X that satisfies

F δ ( X) ≥ n -(⌈∆ -1⌉ -r 0 ) , (21)
we have that

X 0 -X F ≤ nα δ √ ∆ -⌈∆ -1⌉ ,
where α δ = f -1 δ ( 1 n ) . Proof: First, note that due to Assumption 1, f δ (x) takes all the values in ]0, 1[ exactly twice; once with a positive x and once with a negative one. Because of the symmetry, the two have the same modulus; therefore, α δ is well-defined.

Let us denote the singular values of X 0 and X by σ 1 ≥ • • • ≥ σ n and σ1 ≥ • • • ≥ σn , respectively. Define I α as the set of indices i for which σi > α δ . Now, we have that

F δ ( X) = n i=1 f δ (σ i ) = i∈Iα f δ (σ i ) < 1 n <n 1 n =1 + i / ∈Iα f δ (σ i ) ≤1 ≤n-|Iα| < n -|I α | + 1.
On the other hand,

F δ ( X) ≥ n -(⌈∆ -1⌉ -r 0 ); therefore, n -(⌈∆ -1⌉ -r 0 ) < n -|I α | + 1 ⇒ |I α | < (⌈∆ -1⌉ -r 0 ) + 1 ⇒ |I α | ≤ ⌈∆ -1⌉ -r 0 .
This means that at most ⌈∆ -1⌉ -r 0 singular values of X are greater than α δ . Define

H 0 = 0 X 0 X T 0 0 , H = 0 X X T 0 .
In fact, H 0 and H are symmetric matrices that contain the singular values of X 0 and X, respectively, as their n largest eigenvalues and their negatives as the n smallest eigenvalues [START_REF] Horn | Matrix analysis[END_REF]. Next, we apply Weyl's eigenvalue inequality [START_REF] Horn | Matrix analysis[END_REF] as

λ ⌈∆-1⌉+1 (H 0 -H) ≤ λ r0+1 (H 0 ) + λ ⌈∆-1⌉-r0+1 (-H) = σ r0+1 + σ⌈∆-1⌉-r0+1 = σ⌈∆-1⌉-r0+1 ≤ α δ ,
where λ i (•) stands for the i-th largest eigenvalue. This reveals the fact that (X 0 -X) has at most ⌈∆ -1⌉ singular values greater than α δ . Since (X 0 -X) is in the null space of A, Corollary 3 implies that

X 0 -X F ≤ nα δ √ ∆ -⌈∆ -1⌉ .
Corollary 4: For the Gaussian function family given in (8), if [START_REF] Parker | Bilinear generalized approximate message passing[END_REF] holds for a solution X ∈ X , then

X -X 0 F ≤ nδ √ 2 ln n √ ∆ -⌈∆ -1⌉ .
Lemma 4: Let f δ , F δ , X , and, X 0 be as defined in Lemma 3 and assume X δ be the maximizer of F δ (X) on X . Then X δ satisfies [START_REF] Parker | Bilinear generalized approximate message passing[END_REF].

Proof: One can write that

F δ (X δ ) ≥ F δ (X 0 ) ≥ n -r 0 ≥ n -(⌈∆ -1⌉ -r 0 ) .
The first inequality comes from the fact that X δ is the maximizer of the F δ (X), and the second one is true because X 0 has (n -r 0 ) singular values equal to zero; thus, in the summation F δ (X) = n i=1 f δ (σ i ), there are (n -r 0 ) ones. Hence, F δ (X 0 ) ≥ n -r 0 . To see the last inequality, note that 2r 0 < ∆ and ∆ ≤ ⌈∆ -1⌉ + 1. Thus, it can be concluded that 2r 0 < ⌈∆ -1⌉ + 1 which results in 2r 0 ≤ ⌈∆ -1⌉ because r 0 ∈ N. Finally, r 0 ≤ ⌈∆ -1⌉ -r 0 which implies that n -(⌈∆ -1⌉ -r 0 ) ≤ n -r 0 .

Lemma 4 and Corollary 4 together prove that for the Gaussian family,

lim δ→0 argmax{F δ (X) | A(X) = b} = X 0 .
In Theorem 4, we extend this result to all function classes that satisfy Assumption 1.

Theorem 4: Suppose A : R n1×n2 → R m has the ∆spherical property and {f δ } satisfies Assumption 1, and define X , F δ , and X 0 as in Lemma 3. If X δ represents the maximizer of F δ (X) over X , then lim δ→0 X δ = X 0 .

Proof: By combining Lemma 3 and Lemma 4, we obtain that

X 0 -X δ F ≤ nα δ √ ∆ -⌈∆ -1⌉ , (22)
where

α δ = f -1 δ ( 1 n
) . The consequence of Assumption 1 in [START_REF] Cox | Multidimensional Scaling[END_REF] shows that for any ǫ > 0 and 0 < x < 1, one can set δ sufficiently small such that f -1 δ (x) < ǫ. Therefore,

lim δ→0 α δ = lim δ→0 f -1 δ 1 n = 0.
This yields lim δ→0 X 0 -X δ F = 0.

V. NUMERICAL SIMULATIONS

In this section, the performance of the SRF algorithm is evaluated empirically through simulations and is compared to some other algorithms. In the first part of numerical experiments, effects of the algorithm parameters (L, c, and ǫ) in reconstruction accuracy are studied. Next, in the second part, the so called phase transition curve [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF] between perfect recovery and failure is experimentally obtained for the SRF algorithm and is compared to that of the nuclear norm minimization. In the third part of simulations, accuracy and computational load of the SRF algorithm in solving the matrix completion problem are compared to five well known matrix completion algorithms. Finally, in the fourth part, robustness of the SRF against the measurement noise is experimentally verified.

To generate a testing random matrix X ∈ R n1×n2 of rank r, the following procedure is used. We generate two random matrices X L ∈ R n1×r and X R ∈ R r×n2 whose entries are independently and identically drawn from a Gaussian distribution with zero mean and unit variance. Then X is constructed as the product of X L and X R , i.e., X = X L X R . Let A ∈ R m×n1n2 denote the matrix representation of A introduced in (2). In the affine rank minimization problems, all entries of A are drawn independently and identically from a zero-mean, unit-variance Gaussian distribution. Moreover, in the matrix completion simulations, the index set Ω of revealed entries is selected uniformly at random. We denote the result of the SRF algorithm by X and measure its accuracy by SNR rec = 20 log 10 ( X F / X -X F ) in dB, which is referred to as the reconstruction SNR. In addition, by term easy problems, we mean problems in which the ratio m/d r is greater than 3, where d r = r(n 1 + n 2 -r) denotes the number of degrees of freedom in a real-valued rank-r matrix [START_REF] Candés | Exact matrix completion via convex optimization[END_REF]. When this ratio is lower than or equal to 3, it is called a hard problem.

In all experiments, the parameter µ is fixed at 1, and we use a decreasing sequence of δ's according to δ j = cδ j-1 , j ≥ 2, where 0 < c < 1 denotes the rate of decay. The value of δ 1 is set twice as large as the largest singular value of the initial estimate. For the sake of simplicity, square matrices are tested, so n 1 = n 2 = n.

Our simulations are performed in MATLAB 8 environment using an Intel Core i7, 2.6 GHz processor with 8 GB of RAM, under Microsoft Windows 7 operating system.

A. Parameters Effects

Experiment 1. As already discussed in Section III-A, it is not necessary to wait for complete convergence of the internal optimization loop. Instead, a few iterations suffice to only move toward the global maximizer for the current value of δ. Thus, we suggested to do the internal loop for fixed L times. However, the optimal choice of L depends on the aspects of the problem at hand. As a rule of thumb, when the problem becomes harder, i.e., the number of measurements decreases toward the degrees of freedom, larger values of L should be used. Likewise, for easier problems, smaller values of L decrease the computational load of the algorithm, while the accuracy will not degrade very much.

To see the above rule, the affine rank minimization problem defined in (1) is solved using the SRF algorithm, while changing the parameter L. We put n = 30, r = 3, ǫ = 10 -5 , and c = 0.9. The number of measurements changes from 250 to 500 to cover both easy and hard problems. To obtain accurate SNR rec estimates, the trials are repeated 100 times. Fig. 2 shows the effects of changing L from 1 to 10. It can be concluded from Fig. 2 that for easy and hard problems, there is a threshold value for L, which choosing L beyond it can only slightly improves reconstruction SNR. However, in our simulations, we found that increasing the L boosts the computation time almost linearly. For instance, when m = 500 and L = 1, the average computation time is about 0.5 sec, while this time increases to about 1.2 sec for L = 5 and to about 2.2 sec for L = 10.

Experiment 2. The next experiment is devoted to the dependence of the accuracy of the SRF algorithm on the parameter c. In this experiment, the dimensions of the matrix are the same as the in previous experiment, and L and ǫ are fixed to 8 and 10 -5 , respectively. Affine rank minimization and matrix completion problems are solved with two different number of measurements to show the effect on different conditions. c is changed from 0.15 to 0.95, SNR rec 's are averaged on 100 runs. Fig. 3 depicts the reconstruction SNR versus the parameter c for different problems. It is obvious that SNR increases as c approaches 1. However, when c exceeds a critical value, SNR remains almost constant.

Generally, the optimal choice of c depends on the criterion which aimed to be optimized. When accuracy is the key criterion, c should be chosen close to 1, which results in slow decay in the sequence of δ and a higher computational time.

Experiment 3. In this experiment, the effect of ǫ on the accuracy of the algorithm is analyzed. All dimensions and parameters are the same as in the experiment 2 except c and MC, m=550 MC, m=450 RM, m=400 RM, m=300 Fig. 3. Averaged SNRrec of the SRF algorithm as a function of c. Matrix dimensions are fixed to 30 × 30 and r is set to 3. The parameter L and ǫ are set to 8 and 10 -5 , respectively, to have small effect on this analysis. SNR's are averaged over 100 runs. 'MC' and 'RM' denote the matrix completion and affine rank minimization problems, respectively. For two MC problems, m is set to 450 and 550, and for two RM problems, m is set to 300 and 400.

ǫ. c is fixed to 0.9, and ǫ is changed from 10 -1 to 10 -6 . The result of this experiment is shown in Fig. 4. It is seen that after passing a critical value, logarithmic reconstruction SNR increases almost linearly as ǫ decreases linearly in logarithmic scale. Hence, it can be concluded that ǫ controls the closeness of the final solution to the minimum rank solution.

B. Phase Transition Curve

Experiment 4. To the best of our knowledge, the tightest available bound on the number of required samples for the NNM to find the minimum rank solution is two times greater than that of the rank minimization problem [START_REF] Mohan | A simplified approach to recovery conditions for low rank matrices[END_REF]. More precisely, for the given linear operator which has a null space with ∆-spherical section property, (1) has a unique MC, m=550 MC, m=450 RM, m=400 RM, m=300 Fig. 4. Averaged SNRrec of the SRF algorithm as a function of ǫ. Matrix dimensions are fixed to 30 × 30 and r is set to 3. The parameter L and c are set to 8 and 0.9, respectively, to have small effect on this analysis. SNR's are averaged over 100 runs. ǫ is changed from 10 -1 to 10 -6 . 'MC' and 'RM' denote the matrix completion and affine rank minimization problems, respectively. For two MC problems, m is set to 450 and 550, and for two RM problems, m is set to 300 and 400. solution if rank(X 0 ) < ∆/2, while ( 4) and ( 1) share a common solution if rank(X 0 ) < ∆/4. Our main goal, in this experiment, is to show that the SRF algorithm can recover the solution in situations where nuclear norm minimization fails. In other words, this algorithm can get closer to the intrinsic bound in recovery of low-rank matrices. The computational cost of the SRF algorithm will be compared to an efficient implementation of the nuclear norm minimization in the next experiment.

Like compressive sensing literature, the phase transition can be used to indicate the region of perfect recovery and failure [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]. Fig. 5 shows the results of applying the proposed algorithm on the affine rank minimization. A solution is declared to be recovered if reconstruction SNR is greater than 60 dB. The matrix dimension is 40 × 40, ǫ = 10 -5 , L = 6, and c = 0.9. Simulations are repeated 50 times. The gray color of cells indicates the empirical recovery rate. White denotes perfect recovery in all trials, and black shows unsuccessful recovery for all experiments. Furthermore, the thin trace on the figure shows a theoretical bound in recovery of low rank solutions via the nuclear norm minimization found in [START_REF] Oymak | New null space results and recovery thresholds for matrix rank minimization[END_REF]. In [START_REF] Oymak | New null space results and recovery thresholds for matrix rank minimization[END_REF], it is shown that this bound is very consistent to the numerical simulations; thus, we use it for the sake of comparison. One can see in Fig. 5 that there is a very clear gap between this bound and phase transition of the SRF algorithm.

C. Matrix Completion

Experiment 5. The accuracy and computational costs of the proposed algorithm in solving the matrix completion problem are analyzed and compared to five other methods. Among many available approaches, IALM [START_REF] Lin | The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices[END_REF], APG [START_REF] Toh | An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems[END_REF], LMaFit [START_REF] Wen | Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm[END_REF], BiG-AMP [START_REF] Parker | Bilinear generalized approximate message passing[END_REF], and OptSpace [START_REF] Keshavan | Matrix completion from a few entries[END_REF] Gray-scale color of each cell indicates the rate of perfect recovery. White denotes 100% recovery rate, and black denotes 0% recovery rate. A recovery is perfect if the SNRrec is greater than 60 dB. The red trace shows the so called weak bound derived in [START_REF] Oymak | New null space results and recovery thresholds for matrix rank minimization[END_REF] for the number of required measurements for perfect recovery of low rank matrix using the nuclear norm heuristics.

complexity [START_REF] Lin | The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices[END_REF], [START_REF] Toh | An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems[END_REF], while other selected methods are only applicable to the MC setting and exploit other heuristics rather than the nuclear norm to find a low-rank solution.

LMaFit, which is known to be very fast in completing partially observed matrices, uses a nonlinear successive over-relaxation algorithm [START_REF] Wen | Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm[END_REF]. BiG-AMP extends the generalized approximate message passing algorithm in the compressive sensing to the matrix completion and outperforms many state-of-theart algorithms [START_REF] Parker | Bilinear generalized approximate message passing[END_REF]. OptSpace is based on trimming rows and columns of the incomplete matrix followed by truncation of some singular values of the trimmed matrix [START_REF] Keshavan | Matrix completion from a few entries[END_REF]. LMaFit, BiG-AMP, and OptSpace require an accurate estimate of the rank of the solution. MATLAB implementation of OptSpace6 is provided with a function for estimating the rank of the solution, and we use it in running OptSpace. Moreover, LMaFit7 should be initialized with an upperbound on the rank of the solution which, in our numerical experiments, this upperbound is set to 1 2 n. Also, BiG-AMP8 needs a similar upperbound to learn the underlying rank, and we pass 1 2 n as the upperbound to the Big-AMP algorithm too.

IALM 9 , LMaFit, and OptSpace are run by their default parameters except for tol = 10 -9 . For APG 10 , we use default parameters and set tol and mu_scaling to 10 -9 to have the best achieved SNR rec on the same order of other methods. SRF is run with ǫ = 10 -9 , L = 8, and c = 0.95.

Matrix dimensions are fixed to 100 × 100, and r is set to 8, 16, and 32. To see the performance of the aforementioned algorithms, SNR rec and execution time are reported as a function of m/d r for the three values of the rank. Although CPU time is not an accurate measure of the computational costs, we use it as a rough estimate to compare algorithm complexities. Every simulation is run 100 times, and the results are averaged. Fig. 6 demonstrates the results of these comparisons for the three matrix ranks as a function of number of measurements. In comparison to BiG-AMP, while SRF starts completing low-rank matrices with a good accuracy approximately with the same number of measurements when the rank equals to 8, once r increases to 16, it needs smaller number of measurements to successfully recover the solutions. This gap is widen when r = 32. Furthermore, in all simulated cases, SRF has lower running time when compared to BiG-AMP except for starting values of m/d r . SRF also outperforms IALM and APG, which implement NNM, in terms of accuracy, whereas its computational complexity is very close to that of APG. Finally, although the execution time of LMaFit is considerably lower than that of SRF, it needs much larger number of measurements to start recovering low-rank solutions. Note that, here, c is set to 0.95 to accommodate the worst case scenario of hard problems. However, it can be tuned to speed up the SRF method, if the working regime is a priori known.

In summary, the significant advantage of SRF is in solving hard problems where the number of measurements is approaching to d r . Especially, when the matrix rank increases (see Fig. 6(b) and 6(c)), SRF can recover the low-rank solution with at least 20% less number of measurements than other competitors.

D. Noisy Measurements

Experiment 6. Although the SRF algorithm is designed for noiseless measurements, we show experimentally that it is robust against noise. The measurement model can be updated as A(X) = b + ν in the existence of an additive noise vector ν. To check the robustness, we set up an experiment with n = 30, r = 3, c = 0.9, ǫ = 10 -5 . We generate entries of ν as i.i.d. realizations of a zero-mean Gaussian random variable with unit variance. We also vary m from d r to 3d r and average SNR rec over 100 trials. Fig. 7 shows the results of SNR rec vs. the number of measurements for various measurement SNR values (SNR meas = 20 log 10 ( b 2 / ν 2 )). Except for the noiseless case where SNR meas = ∞, we observe that the curves of SNR rec stagnate almost at the same level as SNR meas .

VI. CONCLUSION

In this work, a rank minimization technique based on approximating the rank function and successively improving the quality of the approximation was proposed. We theoretically showed that the proposed iterative method asymptotically achieves the solution to the rank minimization problem, provided that the middle-stage minimizations are exact. We further examined the performance of this method using numerical simulations. The comparisons against five common methods Fig. 6. Comparison of the SRF algorithm with the IALM [START_REF] Lin | The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices[END_REF], APG [START_REF] Toh | An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems[END_REF], LMaFit [START_REF] Wen | Solving a low-rank factorization model for matrix completion by a nonlinear successive over-relaxation algorithm[END_REF], BiG-AMP [START_REF] Parker | Bilinear generalized approximate message passing[END_REF], and OptSpace [START_REF] Keshavan | Matrix completion from a few entries[END_REF] algorithms in terms of accuracy and execution time in completing low-rank matrices. Averaged SNRrec and execution time of all algorithm are plotted as a function of m/dr. Matrix dimensions are fixed to 100 × 100, and r is set to 8, 16, and 32. Trials are repeated 100 times, and results are averaged. No Noise SNR meas =20dB SNR meas =30dB SNR meas =40dB Fig. 7. SNRrec versus number of measurements when there is measurement noise. n = 30, r = 3, c = 0.9, ǫ = 10 -5 , m is changed from dr to 3dr, and results are averaged on 100 runs. reveal superiority of the proposed technique in terms of accuracy, especially when the number of affine measurements decreases toward the unique representation lower-bound. By providing examples, we even demonstrate the existence of scenarios in which the conventional nuclear norm minimization fails to recover the unique low rank matrix associated with the linear constraints, while the proposed method succeeds. lim δ→∞ F δ ( X) = n. Thus, for any ǫ ≥ 0, one can set δ large enough such that F δ ( X) ≥ n-ǫ. Note that for any 1 ≤ i ≤ n, we have that

n -1 + f δ σ i (X δ ) ≥ F δ (X δ ) ≥ F δ ( X) ≥ n -ǫ, or f δ σ i (X δ ) ≥ 1 -ǫ.
This implies that σ i (X δ ) ≤ |f -1 δ (1 -ǫ)| = δ|f -1 (1 -ǫ)|. Hence,

0 ≤ lim δ→∞ σ i (X δ ) δ ≤ f -1 (1 -ǫ) , ∀ 0 < ǫ < 1.
By considering the above inequality for ǫ → 0, we conclude that

lim δ→∞ σ i (X δ ) δ = 0, 1 ≤ i ≤ n.
Using the Taylor expansion, we can rewrite f (•) as f (s) = 1 -γs 2 + g(s),

where γ = -1 2 f ′′ (0) and lim s→0 g(s) s 2 = 0.

In turn, F δ (•) can be rewritten as

F δ (X) = n i=1 f δ σ i (X) = n - γ δ 2 n i=1 σ 2 i (X) + n i=1
g(σ i (X)/δ). [START_REF] Ghasemi | SRF: Matrix completion based on smoothed rank function[END_REF] This helps us rewrite F δ (X δ ) ≥ F δ ( X) in the form

γ δ 2 n i=1 σ 2 i (X δ ) - n i=1 g(σ i (X δ )/δ) ≤ γ δ 2 n i=1 σ 2 i ( X) - n i=1
g(σ i ( X)/δ), or similarly,

σ(X δ ) 2 -σ( X) 2 ≤ n i=1 g σ i (X δ )/δ -g σ i ( X)/δ γ δ -2 ≤ σ(X δ ) 2 γ n i=1 g σ i (X δ )/δ σ i (X δ )/δ 2 + σ( X) 2 γ n i=1 g σ i ( X)/δ σ i ( X)/δ 2 .
Recalling σ(X) 2 = X 2 F , we can write that

X δ 2 F ≤ X 2 F 1 + 1 γ n i=1
g σi( X)/δ σi( X)/δ Application of ( 28) and ( 29) in [START_REF] Candès | Decoding by linear programming[END_REF] results in

lim δ→∞ X δ 2 F ≤ X 2 F . (30)
According to the definition of X, we have X δ 2 F ≥ X 2

F and lim δ→∞ X δ 2 F ≥ X 2 F . Combining this result with [START_REF] Zhang | Theory of compressive sensing via ℓ 1 minimization: A nonrip analysis and extensions[END_REF], we obtain

lim δ→∞ X δ 2 F = X 2 F .
Also, any matrix in null(A) is perpendicular to X since it is the minimum Frobenius norm solution of the A(X) = b. To see this, let A * : R m → R n1×n2 denote the adjoint operator of A and let B : R m → R m denote the inverse of the operator A(A * (•)). Then, similar to the vector case, one can show that X = A * (B(b)) and ∀Z ∈ null(A), X, Z = trace( X T Z) = 0. Thus, X δ 2 F = X 2 F + X δ -X 2 F . In summary, we conclude that lim δ→∞ X δ -X 2 F = 0 which establishes lim δ→∞ X δ = X.

Fig. 1 .

 1 Fig. 1. The SRF Algorithm.

Fig. 2 .

 2 Fig.2. Averaged SNRrec of the SRF algorithm in solving the ARM problem versus L. Matrix dimensions are fixed to 30 × 30, and r is set to 3. The parameter c and ǫ are set to 0.9 and 10 -5 , respectively, to have small effect on this analysis. SNR's are averaged over 100 runs.

Fig. 5 .

 5 Fig. 5. Phase transition of the SRF algorithm in solving the ARM problem. n = 40, ǫ = 10 -5 , L = 6, c = 0.9, and simulations are performed 50 times.Gray-scale color of each cell indicates the rate of perfect recovery. White denotes 100% recovery rate, and black denotes 0% recovery rate. A recovery is perfect if the SNRrec is greater than 60 dB. The red trace shows the so called weak bound derived in[START_REF] Oymak | New null space results and recovery thresholds for matrix rank minimization[END_REF] for the number of required measurements for perfect recovery of low rank matrix using the nuclear norm heuristics.

  r = 16.

  APG(c) r = 32.

2 .

 2 σi(X δ )/δ σi(X δ )/δ

ℓ 0 -norm, not mathematically a vector norm, denotes the number of nonzero elements of a vector.

Note that just r entries of σ(X) are non-zero where r is the rank of the matrix X.

For any finite δ > 0, F δ (•) is not a concave function, and, throughout the paper, a local maximum of F δ (•) denotes a point which is locally and not, at the same time, globally maximum.

To see the definition of subdifferential and subgradient of non-convex functions, refer to[START_REF] Lewis | Convex analysis on the hermitian matrices[END_REF] Sec. 3].

For a convex function, the subdifferential is singleton iff the function is differentiable[START_REF] Rockafellar | Convex Analysis[END_REF].

MATLAB code: web.engr.illinois.edu/ swoh/software/optspace/code.html

MATLAB code: lmafit.blogs.rice.edu/

MATLAB code: sourceforge.net/projects/gampmatlab/

MATLAB code: perception.csl.illinois.edu/matrix-rank/sample code.html

MATLAB code: math.nus.edu.sg/ mattohkc/NNLS.html

ACKNOWLEDGMENT

The authors would like to thank Hooshang Ghasemi for his help in obtaining preliminary results and anonymous reviewers for their helpful comments.

Research Center (ITRC) under contract number 500/11307 and Iran National Science Foundation under contract number 91004600.

APPENDIX A

In this appendix, the closed form least-squares solution of the orthogonal back projection onto the feasible set is derived. Let us cast the affine constraints A(X) = b as A vec(X) = b. The goal is to find the nearest point in the affine set to the result of the j-th iteration, X j . Mathematically,

or equivalently,

where • denotes vector ℓ 2 -norm. By putting y = vec(X)vec(X j ), the problem ( 24) can be easily cast as the following least-squares problem

Let A † = A T (AA T ) -1 be the Moore-Penrose pseudoinverse of A. Then the least-squares solution of (23) will be

, where I denotes the identity matrix, and mat n1,n2 (•) reverses the operation of vectorization, i.e., mat n1,n2 vec(X) = X.

APPENDIX B PROOF OF THEOREM 2

Proof:

To prove lim δ→∞ X δ = X, we first focus on singular values σ i (X δ ). Due to Assumption 1, it is known that Christian Jutten received the Ph.D. and Doctor ès Sciences degrees in signal processing from the Institut National Polytechnique (INP), Grenoble, France, in 1981 and 1987, respectively. In 1982, he became an Associate Professor at the Electrical Engineering Department (1982-1989) of INP, before being Full Professor at University Joseph Fourier of Grenoble, France, in 1989. For 30 years, his research interests have been learning in neural networks, blind source separation, and independent component analysis, including theoretical aspects (separability, source separation in nonlinear mixtures, sparsity, multimodality) and applications (brain imaging, hyperspectral imaging, chemical sensor array, speech). He is author or coauthor of more than 75 papers in international journals, 4 books, 24 keynote plenary talks, and 170 communications in international conferences.

He has been visiting professor at Swiss Federal Polytechnic Institute (Lausanne, Switzerland, in 1989), in Riken labs (Japan, March 1996) and in Campinas University (Brazil, July and August 2010). He was director or deputy director of his lab from 1993 to 2010, especially head of the signal processing department (120 people) and deputy director of GIPSA-lab (300 people) from 2007 to 2010).

He was a scientific advisor for signal and images processing at the French Ministry of Research (1996)(1997)(1998)