A penalized method of alternating projections for weighted low-rank hankel matrix optimization

Weighted low-rank Hankel matrix optimization has long been used to reconstruct contaminated signal or forecast missing values for time series of a wide class. The Method of Alternating Projections (MAP) (i.e., alternatively projecting to a low-rank matrix manifold and the Hankel matrix subspace) is a leading method. Despite its wide use, MAP has long been criticized of lacking convergence and of ignoring the weights used to reflect importance of the observed data. The most of known results are in a local sense. In particular, the latest research shows that MAP may converge at a linear rate provided that the initial point is close enough to a true solution and a transversality condition is satisfied. In this paper, we propose a globalized variant of MAP through a penalty approach. The proposed method inherits the favourable local properties of MAP and has the same computational complexity. Moreover, it is capable of handling a general weight matrix, is globally convergent, and enjoys local linear convergence rate provided that the cutting off singular values are significantly smaller than the kept ones. Furthermore, the new method also applies to complex data. Extensive numerical experiments demonstrate the efficiency of the proposed method against several popular variants of MAP.


Introduction
In this paper, we are mainly interested in the numerical methods for the weighted low-rank Hankel matrix optimization: where H k× is the space of all k × Hankel matrices in the real/complex field with the standard trace inner product, · is the Frobenius norm, M r is the set of matrices whose ranks are not greater than a given rank r , A is given, and W is a given weight matrix W i j ≥ 0. Here (A • B) represents elementwise multiplication (e.g., Hadamard product) between A and B. We note that the size of all matrices involved are of k × . The difficulties in solving (1) are with the low-rank constraint and how to effectively handle a general weight matrix W , the latter of which is often overlooked in existing literature. Our main purpose is to develop a novel, globally convergent algorithm for (1) and its efficiency will be benchmarked against several state-of-the-art algorithms.
In what follows, we first explain an important application of (1) to time series data, which will be tested in our numerical experiment part. We then review the latest advances on algorithms relating to the alternating projection method. We finish this section by explaining our approach and main contributions.

Applications in time series
Problem (1) arises from a large number of applications including signal processing, system identification and finding the greatest common divisor between polynomials [23]. To motivate our investigation on (1), let us consider a complex-valued time series a = (a 1 , a 2 , . . . , a n ) of finite rank [17,Chp. 5]: where P s (t) are a complex polynomial of degree (ν s − 1) (ν s are positive integers) and λ s ∈ C \ {0} are distinct. Define r := ν 1 + . . . + ν m (":=" means "define"). Then it is known [29,Prop. 2.1] that the rank of the Hankel matrix A generated by a: a 1 a 2 · · · a a 2 a 3 · · · a +1 . . . . . . . . . . . . a k a k+1 · · · a n ⎤ ⎥ ⎥ ⎥ ⎦ must be r , where the choice of (k, ) satisfies n = k + − 1 and r ≤ k ≤ n − r + 1. Suppose now that the time series a is contaminated and/or has missing values. To reconstruct a, a natural approach is to computing its nearest time series x by the least squares: where w = (w 1 , . . . , w n ) ≥ 0 are the corresponding weight vector emphasizing the importance of each elements of a. The equivalent reformulation of (3) as (1) is obtained by setting where v is known as the averaging vector of Hankel matrix of size k× (k ≤ ) and √ w is the elementwise square root of w. We note that the widely studied (1) with W i j ≡ 1 corresponds to w i = 1/v i , which is known as the trapezoid weighting. Another popular choice for financial time series is the exponential weights w i = exp(αi) for some α > 0. We refer to [15,Sect. 2.4] for more comments on the choice of weights.
A special type of the time series of (2) arises from the spectral compressed sensing, which has attracted considerable attention lately [4]. In its one dimensional case, a t is often a superposition of a few complex sinusoids: where j = √ −1, r is the model order, ω s is the frequency of each sinusoid, and d s = 0 is the weight of each sinusoid, and τ s ≥ 0 is a damping factor. We note that (4) is a special case of (2) with P s (t) = d s (hence ν s = 1) and λ s = exp(2π jω s − τ s ). If a t is sampled at all integer values from 1 to n, we get a sample vector a ∈ C n . Consequently, the rank of H(a) must be r . However, in practice, only a subset Ω of the sampling points {1, . . . , n} may be observed (possibly contaminated), leading to the question of how to best reconstruct a(t) based on its partial observation a i on Ω. This has led to the Hankel matrix completion/approximation problem of (1), see [4,Sect. II.A] and [2,Sect. 2.1]. A popular choice of W in the spectral compressed sensing is W i j = 1 for all (i, j), resulting in the distance between X and A in (1) being measured by the standard Frobenius norm. In this paper, we assume Assumption 1 W is Hankel and non-negative (i.e., W i j ≥ 0 for all (i, j)).

On alternating projection methods
Low-rank matrix optimization is an active research area. Our short review is only able to focus on a small group of those papers that motivated our research. We note that there are four basic features about the problem (1): (i) X has to be low rank; (ii) X has Hankel structure; (iii) the objective is weighted; and (iv) X may be complex valued. The first feature is the most difficult one to handle because it causes the nonconvexity of the problem. Many algorithms have been developed proposing different ways to handle this low rank constraint. One of the most popular choices is to use the truncated singular value decomposition to project X to be its nearest rank r matrix and we denote the projection by Π M r (X ). This has given rise to the basic Method of Alternating Projections (MAP) (also known as the Cadzow method [1]): Given X 0 ∈ H, update where Π H (·) is the orthogonal projection operator to the Hankel subspace H k× . Despite its popularity in engineering sciences, Cadzow's method can not guarantee the convergence to an optimal solution. Even convergence occurs, not much is known about where it converges to. It has also been criticized for completely ignoring the objective function, see [5,6,8,14]. In particular, the weight matrix W does not enter Cadzow's method at all because the truncated SVD does not admit a closed-form solution under a weighted norm. Gillard and Zhigljavsky [16] proposed to replace Π M r by its diagonally weighted variants and studied how to best approximate the weight matrix W by a diagonal weight matrix. Qi et. al. [25] proposed to use a sequential diagonal weight matrices aiming to get a better approximation to the original weight matrix. Despite the improved numerical convergence, those methods in [16,25] still inherit the essential problem of convergence of Cadzow's method. Recently, Lai and Varghese [20] considered a similar method for a matrix completion problem and established the linear convergence of their method under a kind of "transversality" condition provided that the initial point is close enough to a true rank-r completion. We refer to [7] for a more general transversality condition that ensures a local linear convergence rate of MAP onto nonconvex sets. Alternating projections of Π M r (·) and Π H (·) also play an important role in the class of iterative hard thresholding (IHT) algorithms for spectral compressed sensing. For example, Cai et. al. [3] established the convergence of IHT in the statistical sense (i.e., with high probability) under a coherence assumption on the initial observation matrix A. Although local convergence results (be in the sense of monotonically decreasing [20] or in the statistical sense [3]) may be established for MAP under some conditions, we are not aware of any existing global convergence results mainly due to the nonconvexity of the rank constraint. For the general weighted (1), it appears to be a difficult task to develop a variant of MAP that enjoys both global and local linear convergence properties. We will achieve this through a penalty approach.
Penalty approaches have long been used to develop globally convergent algorithms for problems with rank constraints, see [10,12,13,21,22,27,32,33]. For example, Gao [12] proposed the penalty function p(X ) based on the following observation: where X * is the nuclear norm of X and σ 1 (X ) ≥ · · · σ n (X ) are the singular values of X in nonincreasing order. However, the resulting method, as well as those in [21,22,27], has nothing to do with MAP any more and its implementation is not trivial.

Our approach and main contributions
In this paper, we propose a new penalty function and develop a penalized method whose main step is the alternating projections. We call it the penalized MAP (pMAP). The new penalty function is the Euclidean distance function d M r (X ) from X to M r : Obviously, the original problem (1) is equivalent to We propose to solve the quadratic penalty problem with ρ > 0 being a penalty parameter: By following the standard argument [24,Thm. 17.1] for the classical quadratic penalty method, one can establish that the global solution of (7) converges to that of (1) as ρ approaches infinity provided the convergence happens. However, in practice, it is probably as difficult to find a global solution of (7) as for the original problems. It is hence important to establish the correspondence between the first-order stationary points of (7) and that of (1). This is done in Theorem 1 under a generalized linear independence condition. The remaining task is to efficiently compute a stationary point of (7) for a given ρ > 0. The key observation is that g r (X ) can be represented as the difference of two convex functions, which can be easily majorized (later on its meaning) to get a majorization function g (m) r (X , X ν ) of g r (X ) at the current iterate X ν . We then solve the majorized subproblem: We will show that the update takes the following form: where W (2) := W • W and the division W (2) /(ρ + W (2) ) is taken componentwise. Compared with (5), this update is just a convex combination of the observation matrix A and the MAP iterate in (5). In the special case that W ≡ 0 (which completely ignores the objective in (1)) or ρ = ∞, (9) reduces to MAP. We will analyze the convergence behaviour of pMAP (9). In particular, we will establish the following among others.
(i) The objective function sequence {F(X ν , ρ)} will converge and X ν+1 − X ν converges to 0. Moreover, any limiting point of {X ν } is an approximate KKT point of the original problem (1) provided that the penalty parameter is above certain threshold (see Theorem 2). (ii) If X is an accumulation point of the iterate sequence {X ν }, then the whole sequence converges to X at a linear rate provided that σ r ( X ) σ r +1 ( X ) (see Theorem 3).
Our results in (i) and (ii) provide satisfactory justification of pMAP. It is not only globally convergent, but also enjoys a linear convergence rate under reasonable conditions. Furthermore, we can assess the quality of the solution as an approximate KKT of (1) if we are willing to increase the penalty parameter. Of course, balancing the fidelity term f (X ) and the penalty term is an important issue that is beyond the current paper. The result in (ii) is practically important too. Existing empirical results show that MAP often terminates at a point whose cut-off singular values (σ i , i ≥ r + 1) are significantly smaller than the kept singular values (σ i , i ≤ r ). Such points are often said to have a numerical rank r , but the theoretical rank is higher than r . This is exactly the situation that was addressed in (ii). Those results are stated and proved for real-valued matrices. We will extend them to the complex case, thanks to a technical result (Proposion 2) that the subdifferential of g r (X ) in complex domain can also be computed in a similar fashion as in the real domain. To our best knowledge, this is the first variant of MAP that can handle general weights and enjoys both global convergence and locally linear convergence rate under a reasonable condition (i.e., σ r σ r +1 ). The paper is organized as follows. In the next section, we will first set up our standard notation and establish the convergence result for the quadratic penalty approach (7) when ρ → ∞. Sect. 3 includes our method of pMAP and its convergence results when ρ is fixed. In Sect. 4, we will address the issue of extension to the complex-valued matrices, which arise from (2) and (4). The key concept used in this section is the Wirtinger calculus, which allows us to extend our analysis from the real case to the complex case. We report extensive numerical experiments in Sect. 5 and conclude the paper in Sect. 6.

Quadratic penalty approach
The main purpose of this section is to establish the convergence of the stationary points of the penalized problems (7) to that of the original problem (1) as the penalty parameter ρ goes to ∞. For the simplicity of our analysis, we focus on the real case. We will extend our results to the complex case in Sect. 4. We first introduce the notation used in this paper.

Notation
For a nonnegative matrix such as the weight matrix W , √ W is its componentwise square root matrix ( W i j ). For a given matrix X ∈ C k× , we often use its singular value decomposition (assume k ≤ ) where σ 1 (X ) ≥ · · · ≥ σ n (X ) are the singular values of X and U ∈ C k×k , V ∈ C × are the left and right singular vectors of X . For a given closed subset C ⊂ C k× , we define the set of all projections from X to C by If C is also convex, then P C (X ) is unique. When C = M r , P M r (X ) may have multiple elements. We define a particular element in P M r (X ) that is based on the SVD (10): where U r and V r consist of the first r columns of U and V respectively. Related to the function g r (X ) defined in (6), the function has the following properties by the classical result of Eckart and Young [9]: It follows from [12,Prop. 2.16] that h r (X ) is convex and the subdifferentials of h r (X ) and g r (X ) in the sense of [26,Def. 8.3] are respectively given by where conv(Ω) denotes the convex hull of the set Ω. Finally, we let B (X ) denote the -neighbourhood centred at X .

Convergence of quadratic penalty approach
The classical quadratic penalty methods try to solve a sequence of penalty problems: where the sequence ρ ν > 0 is increasing and goes to ∞. By following the standard argument (e.g., [24,Thm. 17.1]), one can establish that every limit of {X ν } is also a global solution of (1). However, in practice, it is probably as difficult to find a global solution for (13) as for the original problem (1). Therefore, only an approximate solution of (13) is possible. To quantify the approximation, we recall the optimality conditions relating to both the original and penalized problems.
Following the optimality theorem [26,Thm. 8.15], we define the first-order optimality condition of problem (1) and (7): Definition 1 (First-order optimality condition) X ∈ H satisfies the first-order optimality condition of (1) if where λ is the Lagrangian multiplier. Similarly, we say X ν ∈ H satisfies the first-order optimality condition of the penalty problem (7) if We generate X ν ∈ H such that the condition (15) is approximately satisfied: where ν ↓ 0. We can establish the following convergence result.

Theorem 1
We assume the sequence {ρ ν } goes to ∞ and { ν } decreases to 0. Suppose each approximate solution X ν is generated to satisfy (16). Let X be an accumulation point of {X ν } and we assume Then X satisfies the first-order optimality condition (14).
Proof Suppose X is the limiting point of the subsequence {X ν } K . We consider the following two cases. Case 1 There exists an infinite subsequence K 1 of K such that rank(X ν ) ≤ r for ν ∈ K 1 . This would imply ∂ g r (X ν ) = {0}, which with (16) implies P H (∇ f (X ν )) → 0. Hence (14) holds at X with the choice λ = 0.
Case 2 There exists an index ν 0 such that X ν / ∈ M r for all ν 0 ≤ ν ∈ K. In this case, we assume that there exists an infinite subsequence [26, (8.53)]. Therefore, its limit v ∈ ∂d M r ( X ) by the upper semicontinuity. By the assumption (17), we have v / ∈ H ⊥ because v has the unit length. Since H is a subspace, P H (·) is a linear operator. It follows from (16) that Taking limits on {X ν } ν∈K 2 and using the fact Projecting on both sides to H yields Computing the inner product on both sides with P H ((X ν − Π M r (X ν ))/d M r (X ν )), taking limits on the sequence indexed by K 2 , and using the fact P H (ξ ν ) → 0 due to (16), we obtain lim ν∈K 2 We then have Taking limits on both sides of (18) yields which is sufficient for This completes our result. (17) can be interpreted as that any 0 = v ∈ ∂d M r ( X ) is linearly independent of any set of basis of H ⊥ . Therefore, (17) can be seen as a generalization of the linear independence assumption required in the classical quadratic penalty method for a similar convergence result with all the functions involved being assumed continuously differentiable, see [24,Thm. 17.2]. In fact, what we really needed in our proof is that there exists a subsequence {(X ν − Π M r (X ν ))/d M r (X ν )} in Case (ii) such that its limit v does not belong to H ⊥ . That could be much weaker than the sufficient condition (17).

Remark 1 Condition
Theorem 1 establishes the global convergence of quadratic penalty method when the penalty parameter approaches infinity, which drives g r (X ν ) smaller and smaller. In practice, however, we often fix ρ and solve for X ν . We are interested in how far X ν is from being a first-order optimal point of the original problem. For this purpose, we introduce the approximate KKT point, which keeps the first-order optimality condition (15) with an additional requirement that g r (X ) is small enough.
Definition 2 ( -approximate KKT point) Consider the penalty problem (7) and > 0 is given. We say a point X ∈ H is an -approximate KKT point of (1) if

The method of pMAP
This section develops a new algorithm that solves the penalty problem (7), in particular for finding an approximate KKT point of (1). The first part is devoted to the construction of a majorization function for the distance function dist(X , M r ). We then describe pMAP based on the majorization introduced and establish its global and local convergence.

Majorization and DC interpretation
We first recall the essential properties that a majorization function should have. Let θ(·) be a real-valued function defined in a finite-dimensional space X . For a given y ∈ X , we say a function θ (m) The motivation for employing the majorization is that the squared distance function g r (X ) is hard to minimize when coupled with f (X ) under the Hankel matrix structure. It is noted that where the inequality used the fact that Π M r (X ) is a nearest point in M r to X . It is easy to verify that g r (X , Z ) is a majorization function of g r (X ) at Z .
The following way in deriving the majorization is crucial to our convergence analysis. We recall Being the pointwise maximum of linear functions when Z ∈ M r is given, h r (X ) is convex. The convexity of h r (X ) and (12) yields which further implies In other words, g r (X ) can be seen as Difference of Convex functions, the so-called DC function. Using a subgradient is a common way to majorize DC functions, see [12].

The pMAP algorithm
We recall that our main problem is (7). Our first step is to construct a majorized function of F ρ (X ) at the current iterate X ν : (2) . (21) Note that the division is in the sense of componentwise. The subproblem to be solved at iteration ν is where X ν ρ is defined in (21). The last equation in (22) is due to W ρ := ρ + W (2) being Hankel and computing X ν+1 in (22) is equivalent to averaging X ν ρ along its all anti-diagonals. Since A, ρ/(ρ + W (2) ), W (2) /(ρ + W (2) ) are all Hankel matrices (due to Assumption 1), X ν+1 can be calculated through (9).

Algorithm 1 (pMAP)
1: Input data: Matrix A, weight matrix W , penalty parameter ρ, rank r , and the initial X 0 .

Remark 2
Being a direct consequence of employing the majorization technique, the following decreasing property holds: (19)).
A widely used stopping criterion is X ν+1 − X ν ≤ for some small tolerance > 0. We will see below that X ν+1 − X ν approaches zero and hence such convergence check will eventually be satisfied. For our theoretical analysis, we assume that pMAP generates an infinite sequence (e.g., let = 0).

Convergence of pMAP
We have more results on the convergence of pMAP. Theorem 2 Let {X ν } be the sequence generated by pMAP. The following holds.

then X is an -approximate KKT point of (1).
Proof We will use a number of facts to establish (i). The first fact is due to the convexity of f (X ): The second fact is the identity The third fact is due to the convexity of h r (X ) defined in (11) and Π M r (X ) ∈ ∂h r (X ): The last fact is the optimality condition of problem (22): Combining all facts above leads to a sufficient decrease in F ρ (X k ): In the above we also used the fact that X ν+1 − X ν ∈ H. Since the sequence {F ρ (X ν )} is non-increasing and is bounded from below by 0, we have Taking limits on both sides of (26) and using the upper semi-continuity of the projections P M r (·) yields we use the fact that {F ρ (X ν )} is non-increasing to get The first equality holds because g r (X 0 ) = 0 when X 0 ∈ M r . As a result, Therefore, X is an -approximate KKT point of (1).
We note that the first result (i) in Theorem 2 is standard in a majorizationminimization scheme and can be proved in different ways, see, e.g., [32,Thm. 3.7].

Final rank and linear convergence
This part reports two results. One is on the final rank of the output of pMAP and the rank is always bigger than the desired rank r unless A is already an optimal solution of (1). The other is on the conditions that ensure a linear convergence rate of pMAP. For this purpose, we need the following result.
Proposition 1 [11,Thm. 25] Given the integer r > 0 and consider X ∈ IR k× of rank are the singular values of X and u i , v i , i = 1, . . . , r + p are the left and right (normalized) eigenvectors. We assume σ r ( X ) > σ r +1 ( X ) so that the projection operator Π M r (X ) is uniquely defined in a neighbourhood of X . Then Π M r (X ) is differentiable at X and the directional derivative along the direction Y is given by is the tangent subspace of M r at X and Theorem 3 Assume that W > 0 and X be an accumulation point of {X ν }. The following hold.
(i) rank( X ) > r unless A is already the optimal solution of (1).
(ii) Suppose X has rank (r + p) with p > 0. Let σ 1 ≥ σ 2 ≥ · · · ≥ σ k be the singular values of X . Define Under the condition Consequently, the whole sequence {X ν } converges linearly to X.
Proof (i) Suppose X is the limit of the subsequence {X ν } k∈K . We assume rank( X ) ≤ r . It follows from Theorem 2 that Taking limits on both sides of (9) and using the fact that X is Hankel, we get Under the assumption W > 0, we have X = A. Consequently, rank(A) ≤ r , implying that A is the optimal solution of (1). Therefore, we must have rank( X ) > r if the given matrix A is not optimal already.
Moreover, the directional derivative of φ(X ) at X along the direction Y is given by The inequality above holds because Π H (·) is an orthogonal projection to a subspace and its operator norm is 1. The matrices in Proposition 1 have the following bounds.
In the above, we used the fact that ψ(t) := t/(σ r − t) is an increasing function of t for t < σ r . Proposition 1, (29) and (30) imply The second equality above used the fact that the operator norm of Π T Mr ( X ) is not greater than 1 due to T M r ( X ) being a subspace. Since φ(·) is differentiable at X , there exists > 0 such that Therefore, Now we are ready to quantify the error between X ν and X whenever X ν ∈ B ( X ).
Consequently, X ν+1 ∈ B ( X ). Since {X ν } ν∈K converges to X , X ν will eventually falls in B ( X ), which implies that the whole sequence {X ν } will converge to X and eventually converges at a linear rate.

Remark 3 (Implication on MAP)
When the weight matrix W = 0, pMAP reduces to MAP according to (9). Theorem 2(i) implies which obviously implies The decrease property (32) was known in [5, Eq.(4.1)] and was used there to ascertain that MAP is a descent algorithm. Our improved bound (31) says a lightly more that the decrease each step in the function X − Π M r (X ) is strict unless the update becomes unchanged. In this case (W = 0), the penalty parameter is just a scaling factor in the objective, hence the KKT result in Theorem 2(ii) does not apply to MAP. This probably explains why it is difficult to establish similar results for MAP.  [20] is a local result that requires that the initial point of MAP is close enough to M and the rank r assumption of M is also crucial to their analysis, which also motivated our proof. In contrast, our result is a global one and enjoys a linear convergence rate near the limit under a more realistic assumption σ r σ r +1 . One may have noticed that the convergence rate c though strictly less than 1 may be close to 1. This is often numerically observed that MAP often converges slowly. But the more important point here is that in such a situation it ensures that the whole sequence converges. This global convergence justifies the widely used stopping criterion X ν+1 − X ν ≤ .

Extension to complex-valued matrices
The results obtained in the previous sections are for real-valued matrices and they can be extended to complex-valued matrices by employing what is known as the Wirtinger calculus [30]. We note that not all algorithms for Hankel matrix optimization have a straightforward extension from the real case to the complex case, see [6] for comments on some algorithms. We explain our extension below.
Suppose f : C n → IR is a real-valued function in the complex domain. We write z ∈ C n as z = x + jy with x, y ∈ IR n . The conjugatez := x − jy. Then we can write the function f (z) in terms of its real variables x and y. With a slight abuse of notation, we still denote it as f (x, y). In the case where the optimization of f (z) can be equivalently represented as optimization of f in terms of its real variables, the partial derivatives ∂ f (x, y)/∂x and ∂ f (x, y)/∂y would be sufficient. For other cases where algorithms are preferred to be executed in the complex domain, then the Wirtinger calculus [30] is more convenient to use and it is well explained (and derived) in [19]. The R-derivative and the conjugate R-derivative of f in the complex domain are defined respectively by The R-derivatives in the complex domain play the same role as the derivatives in the real domain because the following two first-order expansions are equivalent: Here, we treat the partial derivatives as column vectors and Re(x) is the real part of x.
Note that in the first-order expansion in f (z + Δz) used the conjugate R-derivative. Hence, we define the complex gradient to be ∇ f (z) := 2∂ f /∂z, when it exists. When f is not differentiable, we can extend the subdifferential of f from the real case to the complex case by generalizing (33). In order to extend Theorem 1, we need to characterize ∂d M r (X ) in the complex domain. We may follow the route of [26] to conduct the extension. For example, we may define the regular subgradient of d M r (X ) [26,Def. 8.3] to its complex counterpart by replacing the conjugate-gradient in the first-order expansion in (33) by a regular subgradient. We then define subdifferential through regular subgradients. With this definition in the complex domain, we may extend [26, (8.53)] to derive formulae for ∂d M r (X ). What we needed in the proof of Theorem 1 is (X − Π M r (X ))/d M r (X ) ∈ ∂d M r (X ) when X / ∈ M r . The proof of this result follows a straightforward extension of the corresponding part in [26, (8.53)] and if reproduced here would take up much space. Hence we omit it.
In order to extend the results in Sect. 3, we need the subdifferential of h r (X ) in order to majorize g r (X ). Since h r (X ) is convex, its subdifferential is easy to define. We note that the inner product for the complex matrix space C k× is defined as A, B = Trace (A H B), where A is the Hermitian conjugate, i.e., A H := A T .

Definition 3 The subdifferential ∂h r (X ) is defined as
The following result is really what we needed in order to extend the results in Sect. 3 to the complex domain.

Proposition 2
For any X ∈ C k× , we have P M r (X ) ⊂ ∂h r (X ).
Proof Let Π M r (X ) stand for any element in P M r (X ). It is easy to verify the following identities: We use (11) and (34) to compute This proves the claim.
A direct consequence is that ∂ g r (X ) = X − ∂h r (X ) ⊃ P M r (X ) and the majorization g r (X ) through the subdifferential of h r (X ) holds. The rest of the extension is straightforward and we do not repeat it here.

Numerical experiments
In this section we test two popular problems (time series denoising in real domain and incomplete signal completion in complex domain) to demonstrate the numerical performance of pMAP. The time series denoising problem aims to extract the noiseless data from polluted observations by removing the noise components, while incomplete signal completion problem tries to approximate the missing data in an incomplete complex valued signal. In both numerical experiments, a solver is terminated when any of the following conditions is met here tol, f tol, gtol are set at 10 −5 , 10 −7 and 10 −8 , respectively. A solver will also be terminated if it reaches the maximum iterations, setting at 200. All codes used were written in MATLAB (2019a) and run on a laptop equipped with a 7th Generation Intel Core i5-7200U CPU and 8GB memory card.

Experiment introduction
In the first experiment we compare the proposed pMAP with some leading solvers including Cadzow's method (Cadzow,[14]) and Douglas-Rachford iterations (DRI, [6]) for real-valued time series de-noising. In this test we randomly generate noiseless time series a = (a 1 , a 2 , . . . , a n ) via the following process:  Fig. 1 Plot of F ρ (X ν ) and X ν+1 − X v by pMAP with ρ fixed

Demonstration of convergence
Before coming to the numerical comparison, we first demonstrate the convergence behaviour of Algorithm 1 under different updating strategy of ρ. We plot the sequences of F ρ (X ν ) in Fig. 1a with both W 1 and W 2 . It can be observed that in both cases, the functional value F ρ (X ν ) decreases and converges. We further plot the sequence of X ν+1 − X ν in Fig. 1b using the same example. We find that the sequence of X ν+1 − X ν is also decreasing and converges to zero, which is consistent with Theorem 2.
The behaviours of σ r +1 σ r are shown in Fig. 2 with respect to different ρ updating strategies. In this and later experiments, ρ is initialised as ρ 0 = 10 −2 × m/n 2 where n denotes the length of a time series and m stands for the amount of known observations, which equals to n in this test. As shown in Fig. 2a, σ r +1 σ r approaches zero with increasing ρ, which means g r (X ν ) goes to zero as well. By contrast if ρ is fixed as ρ ν = ρ 0 at each iterate, Fig. 2b shows that σ r +1 σ r decreases much slower than the first strategy. As a result, we will update ρ by ρ ν+1 = 1.1ρ ν at each iterate when ρ ν ≤ n×min(W ), where min(W ) is the minimal weights in W . The behaviour of convergence of Algorithm 1 appears consistence for other choices.

Numerical results
The numerical results are reported in Table 1 including the number of iterations (Iter), cpu time for computation (Time), root of mean square error (RMSE) for each solver which is calculated as wherex = {x 1 , . . . ,x n } is obtained asx = H −1 ( X ) and X is the estimated result from a certain solver. I denotes the index set of all data to be predicted while |I| stands for the size of I. Apparently smaller RMSEs indicate better solution qualities. Because the Cadzow method does not allow users to select arbitrary weights apart from W 1 , we will not report the numerical results for the Cadzow method under W 2 .
We also report success rate (SR) of each solver in Table 1. One instance is successfully de-noised if the relative gap between the RMSE of approximated solution and the best possible RMSE is smaller than a pre-defined threshold, i.e.,

RMSE approx − RMSE best RMSE best ≤ threshold
In this experiment we set threshold = 10%. For any combinations of {n/r /θ }, all data reported in Table 1 are the mean values over 50 randomly generated tests. Our first observation on Table 1 is that when applying W = W 1 , pMAP reports the best results in 3 examples out of 6 while DRI performs the best in the rest 3 examples. In general, Cadzow, DRI and pMAP have very similar performance on estimation accuracy under W 1 because they are MAP-based algorithms.
When the weight matrix is set as W 2 , a significant improvement on the estimation accuracy can be observed for DRI and pMAP comparing to the case W = W 1 . This result matches our expectation because W 2 assumes that all data have equal importance by sharing the same weight, while W 1 implies that data in the middle of a time series are more weighted than the data at both ends.
For all {n/r /θ } combinations, our proposed solver with W 2 always generated the estimation results with lowest RMSEs. It is also important to mention that our pMAP algorithm enjoys the most robust convergence result among all candidate solvers. As a result, we conclude that our proposed pMAP algorithm is competitive and effective in solving real-valued time series denoising problems.

Experiment introduction
In this experiment, we consider the problem of recovering missing values in an incomplete spectral sparse signal. We refer to [2,4] and the references therein for its background in recovering signals which are spectrally sparse via off-grid methodologies. We follow the suggestions in [2] to generate the experiment data a = {a 1 , a 2 , . . . a n } where where j = √ −1, r is the model order, ω s is the frequency of each sinusoid and d s = 0 is the weight of each sinusoid. Both ω s and d s are randomly sampled following uniform distributions, as ω s ∼ U [0, 1) and d s ∼ U [0, 2π). Indexes of missing data are randomly sampled following uniform distribution. In this experiment we introduce three sub-tests focusing on different purposes.
Test a: Incomplete signal recovery without noises In this sub-test we assume only a subset Ω of the sampling points {1, . . . , n} are observed and we aim to recover the signal by estimating the missing data. Here all observed data are noiseless. We use success rate (SR) to measure the performance of the candidate methods in incomplete signal recovery. We say the signal is successfully recovered if wherex is the estimated signal.
In this test, signal length n is set to 499, 999 or 1999, respectively. The percentage of known observations m n is set to 30% or 60% where m stands for the amount of known observations. All combinations of {n/m/r } used in this test can be found in Table 2. To compare the performance of our proposed method, we further introduce three state-of-art solvers in spectral sparse signal recovery including Atomic Norm Minimization (ANM [28]), Fast Iterative Hard Thresholding (FIHT [3]) and Projected Gradient Descent (PGD [2]). W for each solver is defined as if a i+ j−1 is observed, and W i, j = 0 if it is missing.

Test b: Incomplete signal recovery with noises
This sub-test still aims to recover an incomplete spectral sparse signal, but some observed data are noisy while others are noiseless. Here we follow the signal generating process the same as in Test a, however, 1 3 m observations are polluted by random noises, as y i = a i + i where i = θ e t e 2 a . Here e t is a complex stranded normal random variable and the noise level θ is set at 0.2. We assume the index set of polluted observations is known in advance.
The weight matrix W is set as follows. We give polluted observations very small weights (say 1) and noiseless observations are assigned to much larger weights such as 100. The missing data in the signal are given zero weight. This weighting scheme is more reasonable than other choices as we place higher confidence in noiseless data, less confidence in noisy data and no confidence in missing values. It is worth noting that this weight matrix setting is for pMAP only because FIHT and PGD do not support flexible weight choices. The rest settings of this sub-test such as the definition of {n, m, r } are interpreted the same way as in Test a.
Test c: Recover missing data with inaccurately estimated rank In both Test a and Test b, we assume the objective rank r to be a known parameter. However, obtaining the true information of objective rank is quite challenging in many real applications. So in this sub-test we will examine the performances of candidate solvers in recovering the incomplete spectral sparse signals while the objective rank r is incorrectly estimated. Signal length n is set as 3999 and we assume that 30% data in a signal are randomly or accurately observed without noises. True rank r is set as 15 but assumed to be unknown in this test. It means for each solver, we will try different estimated rankr ranges from 6 to 30. Success rates (SR) over 50 instances are reported to measure the performance of each solver.

Test (a)
The numerical results of this test are listed in Table 2 including total iterations (Iter), CPU time in seconds (Time), RMSE and success rate (SR) for each solver. Among all solvers, ANM enjoys the best global convergence result because it is a convex relaxation method. However, the computational cost of ANM is much higher than the rest solvers and it runs out of memories when n is larger than 500. At the same time, it fails to generate better results comparing with other solvers. Hence we do not report its performances in the rest part of this test problem. Although DRI performs slightly better than the Cadzow method in terms of accuracy, both of these two solvers failed to successfully recover any incomplete signals. We further note the fact that in some cases FIHT stopped within a few iteration steps and this behaviour may lead to inferior solutions. Similar behaviours were also reported in another recent research (Fig. 3, [31]). Also with the increasing r , the performance of PGD declines when the ratio between m and r keeps fixed. It is because PGD has some assumptions on the lower bound of m with respect to r (Theorem 2.1, [2]), which may not hold in some cases. On the other hand, pMAP performs the best in 11 cases out of 12 in terms of the SR.
We find that FIHT is more computational efficient than pMAP. This is the result of using the subspace technique in FIHT (Alg. 2 in [3]), which is designed to approximate the matrix projections on to a low rank subspace and therefore reduces the computational time from O(n 3 ) (SVD) to O(nr 2 ). Although this technique can be adapted to our framework so as to reduce the computational cost, there would be no theoretical guarantee for its convergence results. We leave this potential improvement for future researches.
Some applications require to estimate the coefficients of a spectral sparse signal based on the reconstructed signal, including the amplitude d s and the frequency ω s for all s ∈ [1, 2, . . . , r ]. We use the method in [6, Fig.2] to reconstruct the coefficients of recovered signals in this experiment. The reconstruction results are plotted in Fig. 3 over two randomly selected instances. It is easy to observe that when r = 20, FIHT incorrectly estimated most of coefficients with significant errors while both PGD and pMAP could successfully recover most of them. On the other hand when r increases to 40, PGD performed worse than both FIHT and pMAP. There are several coefficients that were accurately estimated by FIHT and pMAP, but failed to be recovered by PGD (those unrecovered data points were shown by pointing arrows in Fig. 3d). Table 3 lists the numerical results for this test including Iter, Time, RMSE and SR for five solvers. Due to the interference of noises, we increase the threshold of success rate to 10 −2 to make sure the numerical results are comparable and meaningful.

Test b)
Experimental results show that our proposed pMAP significantly outperforms the rest four candidate solvers in all tests. All the competitive solvers failed to recover any signals successfully in all cases (i.e., S R = 0 on average). We observed that those solvers were capable of generating points with RMSE at a level of 10 −2 , but encountered extreme difficulty in driving RMSE below 10 −2 . In contrast, the SR of pMAP is at least 0.92. It might be due to the fact that in DRI, FIHT and PGD, the weight of each observations can not be customised, which means these solvers have to give equal weights to both noisy observations and noiseless observations. As a result, their estimation results are significantly affected by noisy observations. This test clearly demonstrates the advantage of pMAP by allowing customized weighting schemes.
Test (c) The numerical results of recovering missing data with inaccurately estimated rank experiment are plotted in Fig. 4 for each solver. Whenr is smaller than 15, success rate for all solvers are zero. It indicates that none of these solvers can recover the incomplete signal successfully when there is a lack of coefficients information. Withr exactly equals to 15, all three methods including FIHT, PGD and pMAP can achieve 100% recovery rate. However whenr > 15, one can expect varying performances of the three solvers. The success rates of both PGD and FIHT gradually decline with the increasingr and they finally reach around 40% whenr = 30. One     the other hand, SR of our proposed pMAP stays at 100% for anyr no smaller than 15, which indicates that pMAP is more robust to the overestimation of objective rank than the other two solvers. One may wonder what has caused the failure of FIHT and PGD in this case. We explain it below. Figure 5 compares the PGD, FIHT and pMAP using a random incomplete signal recovery example in terms of the relative error between the two consecutive iterates x ν+1 and x ν . Figure 5a shows that the failures of both FIHT and PGD in recovering incomplete signal was caused by the non-convergence behaviour. Figure 5b illustrates the distribution of singular values at the final iterates for each algorithm. A noticeable feature is that FIHT terminates too early because it cut off too many singular values that are not negligible, while PGD and pMAP cut off all comparably small eigenvalues. However, PGD suffers from non-monotonic convergence as shown in Fig. 5a. In contrast, the performance of pMAP just fits the situation studied in Theorem 3, which ensures locally linear convergence when the cut-off singular values are negligible comparing to the first 15 largest ones.

Conclusion
In this paper, we studied the problem of approximating a low rank Hankel matrix from the noisy and/or incomplete observation matrix under arbitrary weights. It is very challenging to tackle this problem due to its non-convexity. We introduced the framework of majorization minimization method such that this problem can be tackled iteratively with non-increasing objective function values. We further showed that the subproblem enjoys a closed form solution, which can be efficiently computed. We demonstrated the global optimal convergence property of our approach pMAP by assuming that the penalty parameter goes to infinity. We also showed that this method will at least converge to an -approximate KKT point linearly if the penalty parameter ρ is above a threshold. This method can be extended to tackle complex-valued matrices because the majorization g r (X ) through the subdifferential of h r (X ) holds in the complex-valued case. In the computational experiments for both series denoising and signal completion problems, pMAP usually outperforms other state-of-the-art solvers in terms of approximation accuracy within reasonable computing times.
One of the important topics to be further investigated is whether the computing cost of pMAP can be improved. For example, the subspace optimization technique used in [3] can reduce the computing time significantly compared with partial SVD, which is the major source of computing cost in pMAP. However, this technique at its current form is likely to break the established convergence result of pMAP, simply because it can not guarantee a closed form solution of each subproblem. Another interesting future research topic is extending our proposed pMAP framework to other rank-minimization related problems with similar structures, such as robust matrix completion and robust principal component analysis.
One referee suggested to consider the exact penalty using the distance function itself rather than its squared from, and use the majorization-minimization approach to the penalized problem. While we recognized the benefit of using the squared distance function as well surveyed in [18], we also think the suggested approach is worth serious investigation. The focus of the difficulty is on how to design an efficient majorization function for the distance function. We leave this to our next research topic.
Code availability pMAP v1.0.0 is available under GNU General Public License. The URLs are contained in this published paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.