Lp Quasi-norm Minimization: Algorithm and Applications

Sparsity finds applications in diverse areas such as statistics, machine learning, and signal processing. Computations over sparse structures are less complex compared to their dense counterparts and need less storage. This paper proposes a heuristic method for retrieving sparse approximate solutions of optimization problems via minimizing the ℓp quasi-norm, where 0<p<1. An iterative two-block algorithm for minimizing the ℓp quasi-norm subject to convex constraints is proposed. The proposed algorithm requires solving for the roots of a scalar degree polynomial as opposed to applying a soft thresholding operator in the case of ℓ1 norm minimization. The algorithm’s merit relies on its ability to solve the ℓp quasi-norm minimization subject to any convex constraints set. For the specific case of constraints defined by differentiable functions with Lipschitz continuous gradient, a second, faster algorithm is proposed. Using a proximal gradient step, we mitigate the convex projection step and hence enhance the algorithm’s speed while proving its convergence. We present various applications where the proposed algorithm excels, namely, sparse signal reconstruction, system identification, and matrix completion. The results demonstrate the significant gains obtained by the proposed algorithm compared to other ℓp quasi-norm based methods presented in previous literature.


Motivation
In numerical analysis and scientific computing, a sparse matrix/array is the one with many of its elements being zeros.The number of zeros divided by the total number of elements is called sparsity.Sparse data is often easier to store and process.Hence, techniques for deriving sparse solutions and exploiting them have attracted the attention of many researchers in various engineering fields like machine learning, signal processing, and control theory.
The taxonomy of sparsity can be studied through the Rank Minimization Problem (RMP).It has been lately considered in many engineering applications including control design and system identification.This is because the notions of complexity and system order can be closely related to the matrix rank.The RMP can be formulated as follows: where X ∈ R m×n and M ⊂ R m×n is a convex set.The problem (1) in its generality is NP- hard [1].Therefore, polynomial time algorithms for solving large-scale problems of the form in (1) are not currently known.Hence, recently adopted methods for solving such problems are approximate and structured heuristics.A special case of RMP is the Sparse Vector Recovery (SVR) problem involving ℓ 0 pseudo-norm minimization given by: where x ∈ R n , V ⊂ R n is a closed convex set and �•� 0 counts the number of the non-zero elements of its argument.From the definition of the rank being the number of non-zero singular values of a matrix, it can be easily realized that (1) is a generalized form of (2).
Numerous studies, which will be expounded upon in the subsequent section, have individually addressed effective solution methods for the problems presented in (1) and (2).These approaches utilize Schatten-p and ℓ p quasi-norm relaxations, respectively.However, existing methods in this domain often either assume a predefined structure for the convex set M in (1) or exclusively cater to the specialized case articulated in (2).Consequently, these methods lack comprehensive applicability.Leveraging the inherent relationship between the Schatten-p quasi-norm and the ℓ p quasi-norm of matrix singu- lar values, we endeavor to formulate an efficient heuristic method based on Schatten-p relaxation.This method is devised to address both problems in a unified manner.The proposed approach begins with the introduction of an algorithm for solving the ℓ p quasi- norm relaxation of the SVR problem presented in (2).Subsequently, recognizing that (2) constitutes a specific case of (1), we utilize the developed ℓ p quasi-norm minimization algorithm as a foundational component for constructing the envisaged generalized algorithm for RMPs.

Sparse vector recovery
Given that many signals exhibit sparsity or compressibility, the SVR problem has found widespread applications in fields such as object recognition, classification, and compressed sensing, as evidenced by studies such as [2][3][4].The concept of sparse representation of signals and systems has been extensively discussed in [5], where the authors conducted a comprehensive review of both theoretical and empirical results pertaining to sparse optimization.They also derived the sufficient conditions necessary for ensuring uniqueness, stability, and computational feasibility.Moreover, [5] explores diverse applications of the SVR problem, contending that in certain tasks involving denoising and compression, methods rooted in sparse optimization offer state-of-the-art solutions.
The problem of constructing sparse solutions for undetermined linear systems has garnered significant attention.A survey conducted in [6] comprehensively examined existing algorithms for sparse approximation.The reviewed methods encompassed various approaches, including greedy methods [7,8], techniques rooted in convex relaxation [3, (1) min X∈M Rank(X), (2) 4], those employing non-convex optimization strategies [9,10], and approaches necessitating brute force [11].The authors discussed the computational demands of these algorithms and elucidated their interrelationships.
Sparse optimization problems of the form min f (x) + µg(x) have been extensively explored in the literature, where g(x) serves as a sparsity-inducing function, f represents a loss function capturing measurement errors, and µ > 0 functions as a trade-off param- eter balancing data fidelity and sparsity.In [12], the authors addressed a sparse recovery problem involving a set of corrupted measurements.By defining g(•) as the ℓ 1 norm, they established a sufficient condition for exact sparse signal recovery, specifically the Restricted Isometry Property (RIP).
Motivated by the convergence of the ℓ p quasi-norm to the ℓ 0 pseudo-norm as p → 0 , the problem was extended in [13] by setting g as the ℓ p quasi-norm for p ∈ (0, 1) .The authors presented theoretical results showcasing the ℓ p quasi-norm's capability to recover sparse signals from noisy measurements.Under more relaxed RIP conditions, it was demonstrated that the ℓ p quasi-norm provides superior theoretical guarantees in terms of stability and robustness compared to ℓ 1 minimization.
In [9], the authors considered the problem of SVR via ℓ p quasi-norm minimization from a limited number of linear measurements of the target signal.However, the proposed approach faced limitations due to its higher computational complexity compared to the ℓ 1 norm.In [14], Fourier-based algorithms for convex optimization were leveraged to solve sparse signal reconstruction problems via ℓ p quasi-norm minimization, demon- strating a combination of the construction capabilities of non-convex methods with the speed of convex ones.
An alternative approach for sparse reconstruction was proposed in [15], replacing the non-convex function with a quadratic convex one.Furthermore, [16] introduced an Alternating Direction Method of Multipliers (ADMM) [17] based algorithm enforcing both sparsity and group sparsity using non-convex regularization.Additionally, [18] proposed an iterative half-thresholding algorithm for expedited solutions of ℓ 0.5 regulariza- tion.The authors not only established the existence of the resolvent of the gradient of the ℓ 0.5 quasi-norm but also derived its analytic expression and provided a threshold- ing representation for the solutions.The convergence of this iterative half-thresholding algorithm was studied in [19], demonstrating its convergence to a local minimizer of the regularized problem with a linear convergence rate.
Conditions for the convergence of an ADMM algorithm aimed at minimizing the sum of a smooth function with a bounded Hessian and a non-smooth function are established in [20].In [21], the convergence of ADMM is analyzed for the minimization of a non-convex and potentially non-smooth objective function subject to equality constraints.The derived convergence guarantee extends to various non-convex objectives, encompassing piece-wise linear functions, ℓ p quasi-norm, and Schatten-p quasi-norm ( 0 < p < 1 ), while accommodating non-convex constraints.Several works have explored the ℓ 1−2 relaxation objective, defined the difference between ℓ 1 and ℓ 2 norms, i.e., ℓ 1 − ℓ 2 , with [22] providing a theoretical analysis on SVR through weighted ℓ 1−2 mini- mization when partial support information is available.Recovery conditions for exact SVR within a ℓ 1−2 objective framework are derived in [23,24], along with references therein, establishing the theoretical foundation for ensuring accurate SVR outcomes.

Rank minimization
In [25], the authors sought to determine the least order dynamic output feedback, utilizing the formulation akin to (1), capable of stabilizing a linear time-invariant system.Their approach involved minimizing the trace, as opposed to the rank, resulting in a Semi-Definite Program (SDP) amenable to efficient solution techniques.Notably, their solution was specifically applicable to symmetric and square matrices.Building upon this work, [26] introduced a generalization of the aforementioned approach.This extension involved replacing the rank in the objective function with the summation of the singular values of the matrix, commonly known as the nuclear norm.The authors demonstrated that this modification yields the convex envelope of the non-convex rank objective, reducing to the original trace heuristic when the decision matrix assumes the form of a symmetric Positive Semi-Definite (PSD) matrix.
In [27], an alternative heuristic based on the logarithm of the determinant was introduced as a surrogate for rank minimization within the subspace of PSD matrices.The authors demonstrated that this formulation could be effectively solved through a sequence of trace minimization problems.In a related study, [28] delved into existing trace and log determinant heuristics, exploring their applications for computing a lowrank approximation in various scenarios.Specifically, the applications encompassed obtaining simple data models with interpretability by approximating covariance matrices for a given dataset.
Drawing inspiration from the success of the ℓ p quasi-norm (0 < p < 1) for sparse sig- nal reconstruction, an alternative method aims to enforce low-rank structure using the Schatten-p quasi-norm.This norm is defined as the ℓ p quasi-norm of the singular val- ues.In [29], the authors addressed the matrix completion problem, which involves constructing a low-rank matrix based on a subset of its entries.Instead of minimizing the nuclear norm, they proposed a Schatten-p quasi-norm formulation and investigated its convergence properties.To enhance the robustness of the solution, [30] combined the Schatten-p quasi-norm for low-rank recovery with the ℓ p quasi-norm (0 < p ≤ 1) of pre- diction errors on the observed entries.The authors introduced an algorithm based on ADMM, which demonstrated superior numerical performance compared to other completion methods.In a non-convex approach for matrix optimization problems involving sparsity, [31] developed a technique using a generalized shrinkage operation.This method enhances the separation of moving objects from the stationary background by decomposing video into low-rank and sparse components, presenting advantages over the convex case.

Contributions
In spite of the commendable performance exhibited by the array of algorithms outlined in Sects.1.2.1 and 1.2.2, each designed to address different relaxations of (1) and (2), it is essential to acknowledge their problem-specific nature, primarily grounded in the specific structural attributes of the convex constraint sets they address.This issue of specialization results in a lack of generality across problem domains.
In this paper, we present a versatile algorithm grounded in the principles of projections onto constraint sets.A distinctive feature of this approach lies in its minimal reliance on problem-specific structural constraints, prioritizing the foundational characteristic of closed convexity.The works [32,33] delve into a comprehensive exploration, analyzing the intrinsic attributes of the projection operation onto constraint sets.While the former addresses the issue without incorporating a crucial coupling condition for polynomial equations, the latter assumes prior knowledge of the projection technique for each given point on ℓ p balls.
Initially, we propose an ADMM based algorithm, termed as ℓ p Quasi-Norm ADMM (pQN-ADMM), designed to solve the ℓ p quasi-norm relaxation of (2).At each itera- tion, the pivotal operation involves computing Euclidean projections onto specific convex and non-convex sets.Notably, the algorithm exhibits two key properties: 1) Its computational complexity aligns with that of ℓ 1 minimization algorithms, with the additional task of solving for the roots of a polynomial; 2) It does not necessitate a specific structure for the convex set.
Subsequently, we extend the application of the proposed algorithm to address the relaxation of (1) by embracing the Schatten-p quasi-norm.In this extension, we leverage the equivalence between minimizing the ℓ p quasi-norm of the vector of singular values and minimizing the Schatten-p quasi-norm.Our study encompasses the following numerical instances: 1 An example employing SVR, wherein the primary objective is the recovery of the sparsest feasible vector from given realizations.2 A matrix completion example, where the overarching goal is the reconstruction of an unknown low-rank matrix based on a limited subset of observed entries.3 Addressing a time-domain system identification problem, specifically tailored for minimum-order system detection.
Our numerical results compellingly showcase the competitiveness of pQN-ADMM when bench-marked against several state-of-the-art baseline methods.
Conclusively, given the inherent reliance of the derived algorithm on a convex projection step in each iteration, our endeavor is directed towards the formulation of an expedited algorithm accompanied by a rigorous mathematical convergence guarantee.Focusing on a subset of problems where the constraint set manifests as a polytope, we leverage principles from the Proximal Gradient (PG) method to formulate a rapid algorithm.The convergence of this algorithm is established with a rate of O( 1K ) , where K denotes the iteration budget assigned to the algorithm.

Notation
Unless otherwise specified, we denote vectors with lowercase boldface letters, i.e., x , with i-th entry as x i , while matrices are in uppercase, i.e.X , with (i, j)-th entry as x i,j .For an integer n ∈ Z + , [n] � ={1, . . ., n} . 1 represents a vector of all entries equal to 1, while 1 G (.) is an indicator function to the set G , i.e., it evaluates to zero if its argument belongs to the set G and is +∞ otherwise.
For a vector x ∈ R n , the general ℓ p norm is defined as: where, we let x be the well-known Euclidean norm, i.e., p = 2 .When 0 < p < 1 , the expression in ( 3) is termed as a quasi-norm satisfying the same axioms of the norm except the triangular inequality making it a non-convex function.
For a matrix X , X represents the spectral norm, which is defined as the square root of the maximum eigenvalue of the matrix X H X .X H refers to the complex conjugate transpose of X , denoted as X ⊤ .On the other hand, .f signifies the Frobenius norm of a matrix.
The Schatten-p quasi-norm of a matrix X is defined as: where σ i (X) is the i-th singular value of the matrix X .We utilize the * subscript in (4) to differentiate the matrix Schatten-p quasi-norm from vector ℓ p case defined in (3).When p = 1 , (4) yields the nuclear norm which is the convex envelope of the rank function.
Throughout the paper, we consider a non-convex relaxation for the rank function, specifically p = 1/2.We define the ceiling operator, denoted as ⌈•⌉ , the vectorization operator vec(X) ∈ R mn , representing the vector obtained by stacking the columns of the matrix X ∈ R m×n , and the Hankel operator Hankel(.), producing a Hankel matrix from the provided vector arguments.We define the sign operator, denoted as sign(•) , which outputs -1, 0, or 1 cor- responding to a negative, zero, or positive argument, respectively.

Problem formulation
This section develops a method for approximating the solution of (2) using the following relaxation: where p ∈ (0, 1] and V is a closed convex set.Problem ( 5) is convex for p ≥ 1 ; hence, can be solved to optimality efficiently.However, the problem becomes non-convex when p < 1 .We present a gradient-based algorithm and consequently, it may not always con- verge to a global optimum solution but only to a stationary point.An epigraph equivalent formulation of ( 5) is obtained by introducing the variable t Let X ⊂ R 2 denote the epigraph of the scalar function |x| p , i.e., X = {(x, t) ∈ R 2 : t ≥ |x| p } , which is a non-convex set for p < 1 .Then, (6) can be cast as: (3) ADMM, as introduced in [17], leverages the inherent problem structure to partition the optimization process into simpler sub-problems, which are solved iteratively.To achieve this, auxiliary variables y = [y i ] i∈[n] and z = [z i ] i∈[n] are introduced, leading to an ADMM reformulation of the problem defined in (7): The dual variables associated with the constraints x = y and t = z are and θ , respec- tively.Throughout the paper, the colons in the constraints of an optimization problem serve as a means to associate the constraint (appearing on the left side of the colon) with its corresponding Lagrange multiplier (found on the right side of the colon).The Lagrangian function corresponding to (8) augmented with a quadratic penalty on the violation of the equality constraints with penalty parameter ρ > 0 , is given by: Considering the two block variables (x, t) and (y, z) , ADMM consists of the following iterations: Given the augmented Lagrangian function expressed in (9), it is evident from (10) that the variables x and t are iteratively updated by solving the following non-convex problem: Exploiting the separable structure of ( 14), one immediately concludes that ( 14) can be split into n independent 2-dimensional problems that can be solved in parallel, i.e., for each i ∈ [n]: (7) where � X (.) denotes the Euclidean projection operator onto the set X .Furthermore, ( 9) and (11) imply that y and z are independently updated as follows: Algorithm 1 summarizes the proposed ADMM algorithm.It is clear that z , , and θ merit closed-form updates.However, updating (x, t) requires solving n non-convex problems.Our strategy for dealing with this issue is presented in the following section.

Non-convex projection
In this section, we present the method used to tackle the non-convex projection problem required to update x and t.
As it is clear from ( 15), x and t can be updated element-wise via performing a projec- tion operation onto the non-convex set X , one for each i ∈ [n] .The n projection prob- lems can be run independently in parallel.We now outline the proposed idea for solving one such projection, i.e., we suppress the dependence on the index of the entry of x and t .For (x, t) ∈ R 2 , � X (x, t) entails solving: If t ≥ |x| p , then trivially � X (x, t) = (x, t) .Thus, we focus on the case in which t < |x| p .The following theorem states the necessary optimality conditions for (18).(15) Theorem 1 Let t < |x| p , and (x * , t * ) be an optimal solution of (18).Then, the following properties are satisfied: Proof We prove the statements by contradiction as follows: (a) Suppose that sign(x * ) � = sign(x) , then: i.e., (x * − x) 2 >(0− x) 2 .Hence, g(x * , t * ) − g(0, t * )>0 .Moreover, the feasibility of (x * , t * ) implies that t * > 0 .Thus, (0, t * ) is feasible and attains a lower objective value than that attained by (x * , t * ) .This contradicts the optimality of (x * , t * ).(b) Assume that t * < t .Then: Furthermore, by the feasibility of (x * , t * ) , we have |x * | p ≤ t * < t .Thus, (x * , t) is feasible and attains a lower objective value than that attained by (x * , t * ) .This con- tradicts the optimality of (x * , t * ).(c) Suppose that |x * | p < t , i.e., We now consider two cases, x > 0 and x < 0 .First, let x > 0 .Then, we have by (a) and (21) that 0 < x * < t 1 p .Since t < |x| p , i.e., (x, t) / ∈ X , therefore t 1 p < x and hence, 0 < x * < t 1 p < x .Pick x 0 > 0 such that |x 0 | p = t , i.e., x 0 = t 1 p .Then clearly, x * < x 0 < x .Thus, we have: where the last inequality follows the just proven identity that x * < x 0 < x .Moreo- ver, we have |x 0 | p = t ≤ t * by (b) .Thus, (x 0 , t * ) is feasible and attains a lower objec- tive value than that attained by (x * , t * ) .This contradicts the optimality of (x * , t * ) .
( Furthermore, the feasibility of (x * , t 0 ) follows trivially from the choice of t 0 .Thus, (x * , t 0 ) is feasible and attains a lower objective value than that attained by (x * , t * ) .This contradicts the optimality of (x * , t * ).This concludes the proof.
We now make use of the fact that for (18), an optimal solution (x * , t * ) satisfies t * = |x * | p and hence, (18) reduces to solving: The first order necessary optimality condition for (24) implies the following: By the symmetry of the function |x| p , without loss of generality, assume that x * > 0 and let 0 < p = s q < 1 for some s, q ∈ Z + .A change of variables a q = x * plugged in (25) shows that finding an optimal solution for (18) reduces to finding a root of the following scalar degree 2q polynomial: To determine � X (x, t) , the objective is to find a root denoted as a * for the polynomial in (26), while ensuring that the pair (a * q , a * s ) minimizes the function g(x, t).Algorithm 2 provides a summary of the method employed to address problem (18).

Convex projection
The convex projection for y-update in ( 16) can be formulated as the following convex optimization problem: Convex problems can be solved by a variety of contemporary methods including bundle methods [34], sub-gradient projection [35], interior point methods [36], and ellipsoid methods [37].The efficiency of optimization techniques relies mainly on exploiting the structure of the constraint set.As discussed in Sect.1.3, our objective is to address the problem outlined in (5) with minimal assumptions on the set V .Our only requirement is that V is a closed and convex set.Nevertheless, if feasible, one should capitalize on the inherent structure of V to potentially streamline the computational complexity involved in solving (27).

Rank minimization algorithm
We consider the problem in (1) and propose a method for approximating its solution efficiently.The Schatten-p heuristic of (1) can be written as: where L = min{m, n} and σ i (X) is the ith singular value of X .In the scenario where p = 1 , (28) represents a convex problem, akin to the nuclear norm heuristic.We now consider a non-convex relaxation, specifically for the case where 0 < p < 1 .The problem in (28) attains an epi-graph form: . Defining the epi-graph set Y for the function σ (X) , where (27) where , θ are the dual variables associated with X and t respectively.Similar to (9), the Lagrangian function associated with (31), augmented with a quadratic penalty for the equality constraint violation with a parameter ρ > 0 , can be represented as: where Tr{.} is the trace operator.Given the 2-tuples (X, t) and (Y, z ), the ADMM itera- tions are as follows:

(X, t) Update
By completing the square and employing some straightforward algebraic manipulations, it can be demonstrated that the problem described in (33) is equivalent to: where For simplicity, we will omit the iteration index k.Let's assume that X = P Q ⊤ and X = U V ⊤ represent the Singular Value Decomposi- tion (SVD) of X and X , respectively.Here, and are diagonal matrices with the singu- lar values associated with X and X , while P , U , Q , and V are unitary matrices.Following the steps in [38, Theorem 3], we can express the first term of (38) as: (32) where (a) is because with I L×L being an identity matrix of size L , and exploiting the circular property of the trace while (b) holds is from the main result of [39].In order to make X − Xk 2 f achieve its derived lower bound, we set P = U and Q = V.
The problem in (38) is then equivalent to: are the vectors of singular values of the matri- ces X and X respectively.The optimal solution X * for ( 38) can be determined by first finding the optimal x * for (40), and then obtaining X * = U * V ⊤ , where * = diag(x * ) and diag(.)denotes an operator that transforms a vector into its corresponding diago- nal matrix.Given that the problem in ( 40) is separable, we will proceed by omitting the index i and focus solely on solving: It can be realized that ( 41) is similar to (18), hence, its optimal solution can be found by applying Algorithm 2.

(Y, z) update
Upon updating (X, t) with and θ held constant, the problem in ( 34) can be reformu- lated as: which is clearly a convex optimization problem representing the projection of the point X k+1 + k ρ on the set M and can be solved by various known class of algorithms as dis- cussed in Sect.3.3.Following the update of Y , the update for z in (35) is as follows: which results in a closed-form solution for z k+1 = t k+1 + θ k −1 ρ .

Proximal gradient algorithm
The pQN-ADMM algorithm adeptly handles the ℓ p relaxation of (2), refraining from assuming any specific structure for V beyond its closed and convex nature.Primarily, the algorithm hinges on the computation of Euclidean projections onto V , as outlined in (27).
In this section, we consider a sub-class of problems with a specific structure for the convex set of the form V = {x : f (x) ≤ 0} , where f (x) is a convex function with Lip- schitz continuous gradient.i.e., f is L-smooth: ∇f (x) − ∇f (y) ≤ L x − y for all x, y ∈ R n .Specifically, in order to solve: (40) min we aim to develop an efficient algorithm with some convergence guarantees for the following Lagrangian relaxation: where µ ≥ 0 is the dual multiplier that captures the trade-off between solution sparsity and fidelity.It is imperative to acknowledge that ( 44) and ( 45) exhibit a relationship, albeit not being strictly equivalent.
A canonical problem for the regularized risk minimization has the following form: where h is an L-smooth loss function, and g represents the regularizer term.In cases where both g and h exhibit convexity, the Proximal Gradient (PG) algorithm [40] can iteratively compute a solution to (46) through PG steps.
where prox g/ (.) , for some constant .When g is convex, the proximal map prox g/ is well-defined, thus, the PG step can be computed.
In comparing both ( 45) and ( 46), it is observed that the convexity assumption of in (46) is not met for x p p in (45).When the regularizer is a continuous non-convex function, the proximal map prox g/ may not exist, and computing it in closed form becomes a challenging task.
On the contrary, in the case of x p p , leveraging similar reasoning as employed in the non-convex projection step introduced in Sect.3.2, our objective is to derive an analytical solution that can be efficiently computed.Specifically, assuming p ∈ (0, 1) is a positive rational number, the proposed method for computing the proximal map of x p p involves finding the roots of a polynomial of order 2q, where q ∈ Z + such that p = s/q for some s ∈ Z + .
Since f is L-smooth, for all x, y ∈ R n , we have: Given x k , replacing f (x) with the upper bound in (48) for y = x k , the prox-gradient operation naturally arises as follows: (44) min x �x� p p , s.t.f (x) ≤ 0, (45) By completing the square, (49) yields to: 50) can be rewritten as: which is clearly a separable structure in the entries of x .Therefore, for each i ∈ [n] , we have: where ḡ : R → R + such that ḡ(t) = |t| p for some positive rational p ∈ (0, 1).
Next, we consider a generic form of ( 52), i.e., given some t ∈ R , we would like to compute: The first-order optimality condition ( 53) can be written as: Using similar arguments as in Sect.3.2, we can conclude that the optimal solution t * attains the property that sign(t * ) = sign( t) .Without loss of generality, exploiting the symmetry of the function ḡ , we only consider the case when t > 0 ; hence, the optimal solution t * is the smallest positive root of the following polynomial: Similar to (26), suppose 0 < p = s q < 1 for some positive integers s and q.By employing the variable transformation a (t * ) 1 q , the optimality condition in ( 55) is simplified to the task of finding the roots of a polynomial of degree 2q: (50) x k+1 = arg min (56) a 2q − ta q + 2s qµL a s = 0.
In to solve (44) effectively, we will employ Algorithm 3, which implements the non-convex inexact Accelerated Proximal Gradient (APG) descent method as presented in [41, Algorithm 2].In summary, Algorithm 3 is designed to tackle composite problems of the form in (46), making the assumptions that h is L-smooth and g is a proper lower-semicontinuous function such that F h + g is bounded from below and coer- cive.This means that lim �x�→∞ F (x) = +∞ .It is important to note that neither h nor g are required to be convex.Algorithm 3 can be summarized as follows: • An extrapolation y k is generated as introduced in [42] for the APG algorithm (step 3).• Steps 4 through 9 encompass a mechanism for a non-monotone update of the objective function.Specifically, F (y k ) undergoes scrutiny concerning its relation to the maximum among the most recent l objective values.Step 9 is responsible for adjusting the gradient step accordingly.This adjustment occasionally allows y k to increase the objective, resulting in a situation where F (y k ) becomes lower than the maximum objective value observed in the latest l iterations.• Steps 11 and 12 represent the solution of the PG step using the non-convex projection method.
In the next part, we show that Algorithm 3 converges to a critical point and it exhibits a convergence rate of O( 1 K ) , where K is the iteration budget that is given to the algorithm.

Numerical results
In this section, we present numerical examples to illustrate the application of the pQN-ADMM algorithm, as expounded in Algorithm 1, and the non-convex projection method delineated in Algorithm 2. Within each of the ensuing examples, we conduct comparative analyses with the convex ℓ 1 relaxation solution, achieved through the use of the MOSEK solver [50], and alternative ℓ 0.5 -based solutions previously proposed in the literature.
The degree of the polynomial for which the roots are determined during the non-convex projection step depends on the value of q in the context of p = s q .It might lead one to speculate that the computational complexity of the non-convex projection step is contingent on the specific value of p, suggesting that lower values of p result in slower algorithm performance.In order to explore this aspect, we systematically performed the non-convex projection step 200 times on a vector of 1024 elements, as part of a sparse vector reconstruction example.Throughout this process, we systematically varied the values of the parameter p, considering a range of p values, specifically p ∈ { 1 2 , 1 3 , 1 4 , 1 5 } .The average time to perform the non-convex projection for the entire vector, where the roots of ( 26) for each p are computed using the "root" command in MATLAB, is observed to be nearly constant, approximately 0.03 s.Furthermore, our numerical experiments in this particular example indicated that for p ∈ { 1  3 , 1 4 , 1 5 } , no substantial improve- ment over the ℓ 0.5 case was observed.As a result, these cases are currently undergoing further investigation and are not included in the numerical results section.

Sparse vector recovery (SVR)
In this section, we implement a sparse vector reconstruction problem and compare the solution of the pQN-ADMM algorithm with the ℓ 1 convex relaxation along with an ℓ 0.5 relaxation solution and Linear Approximation for Index Tracking (LAIT), as presented in [51] and [52], respectively.
Let n = 2 10 and m = n/4 , randomly construct the sparse binary matrix, M ∈ R m× n 2 , with a few number of ones in each column.The number of ones in each column of M is generated independently and randomly in the range of integers between 10 and 20, and their locations are randomly chosen independently for each column.Let U = [M, −M] , which is the vertical concatenation of the matrix M and its negative.Following the same setup in [53], the column orthogonality in U is not satisfied.Let x opt ∈ R n be a reference signal with �x opt � 0 = ⌈0.2n⌉, where the non-zero locations are chosen uniformly at random with the values following a zero mean, unit variance Gaussian distribution.Let v = Ux opt + n be the allowable measurement, where n ∈ R m is a Gaussian random vector with zero mean and co-variance matrix σ 2 I m×m , where I is the identity matrix.The sparse vector is reconstructed from v by solving (5) with V = {x : �Ux − v�/�v� − ǫ ≤ 0} , where ǫ = 3σ �v� .All the algorithms are terminated if �x k − x k−1 �/�x k−1 � ≤ 10 −4 or a budget of 200 iterations is consumed.
Figure 1 depicts the correlation between sparsity levels and noise variances concerning solutions derived through ℓ 1 norm minimization, ℓ 0.5 , LAIT, and pQN- ADMM techniques.A threshold of 10 −6 was imposed, designating entries of the solution vector as zero if they fell below this threshold.The reported outcomes are based on the average results obtained from 20 independently conducted random iterations.Notably, it becomes evident that the pQN-ADMM algorithm consistently yields solutions with higher sparsity levels in comparison to its counterpart baseline methods, across a range of σ 2 values.As σ 2 increases, the sparsity level for all approaches decreases, attributable to the heightened scarcity of information pertaining to the original signal within the realization vector, thereby compromising the precision of the reconstruction process.

Rank minimization problem (RMP)
Within this section, our primary focus is directed towards the exploration of the pQN-ADMM algorithm within the RMP framework, as presented in Sect. 4. We commence by engaging in a matrix completion scenario, presenting an extensive comparative analysis pitting the pQN-ADMM algorithm against various baseline methods rooted in the Schatten-p quasi-norm framework.
Additionally, we delve into a time domain system identification example.Notably, we restrict our comparative analysis to the convex nuclear norm.This singularity in focus arises from the unique constraint nature of the problem at hand, specifically the Hankel constraint.To the best of our knowledge, there are no other Schatten-p-based algorithms capable of addressing constraints of this specific nature in the proposed formulation.This serves to underscore the remarkable versatility of the pQN-ADMM algorithm in handling a broad spectrum of constraints, be they within the vector or matrix domain.

Matrix completion
In this section, we apply our algorithm (pQN-ADMM) to a matrix completion example and compare the result to the Matrix Iterative Re-weighted Least Squares (Matrix-IRLS) [54,55], truncated Iterative Re-weighted unconstrained Lq (tIRucLq) [56] and Iterative Re-weighted Least Squares (sIRLS-p & IRLS-p) [57] algorithms.The matrix completion problem is a special case of the low-rank minimization where a linear transform takes a few random entries of an ambiguous matrix X ∈ R m×n .Given only Fig. 1 The influence of noise variance on the sparsity of solutions generated by ℓ 1 norm, ℓ 0.5 , LAIT, and the pQN-ADMM these entries, the goal is to approximate X and find the missing ones.The matrix com- pletion problem with low-rank recovery can be approximated by, where A : R m×n → R q is a linear map with q ≪ mn and b ∈ R q .To facilitate the appli- cation of the aforementioned algorithms, the linear transform A(X) will be reformulated as Avec(X) , where A ∈ R q×mn and vec(X) ∈ R mn represents a vector obtained by stack- ing the columns of the matrix X.
A random matrix M ∈ R m×n with rank r is created using the following method: 1) 2) The entries of both M L and M R are i.i.d Gaussian random variables with zero mean and unit variance.Let M = M + Z , where Z ∈ R m×n is a Gaussian noise with each entry being an i.i.d Gaussian random variable with zero mean and variance σ 2 .The vector b is then created by selecting random q elements from vec( O M) .Since b = Avec( O M) , one can easily construct the matrix A which is a sparse matrix where each row is composed of a value 1 at the index of the corresponding selected entry in the vector b while the rest are zeros.We set m = n = 100 , r = 5 and p = 0.5 .Let d r = r(m + n − r) denotes the dimension of the set of rank r matrices and define s = q mn as the sampling ratio.We assume that s = 0.195 which yields to q = 1950 .It can be realized that d r q < 1 .We set σ = 0.1 , ǫ = 10 −3 , and let the algorithms terminate if a budget of 1000 iterations is reached.To compare the solutions across different algorithms, where X * represents the solution for ( 59 In Fig. 2a, b, we report the average RFD and REtS values for all the algorithms.Despite that, all the baselines are designed to exploit the specific structure of the matrix completion problem, described in (59), while the proposed pQN-ADMM doesn't, it is competitive against them all.This in turn shows the effectiveness of the pQN-ADMM algorithm in solving the rank minimization problems without requiring any prior information about the structure of the associated convex set.

Time domain system identification
We consider a stable Single Input Single Output (SISO) system operating in discrete time, wherein the input vector u ∈ R T corresponds to a temporal span denoted by T, represent- ing the number of input samples.The system is characterized by an impulse response consisting of a fixed number of samples denoted as n.The resultant output of the system is represented by y ∈ R m .However, in practical scenarios, only noisy realizations, denoted as ŷ , are observable.This realization is expressed as ŷ =y + z = h ⊛ u + z , where h ∈ R n sig- nifies the system's original impulse response, z ∈ R m is a random vector with entries drawn independently from a uniform distribution within the range [−0.25, 0.25] , and ⊛ denotes the convolution operator.
Exploiting the window property of convolution, which asserts that m = n + T − 1 , we establish the relationship among the components u i , h i , and y i through the linear convolu- tion relation y i = ∞ j=−∞ h j u i−j .Herein, u i , h i , and y i represent the ith components of the vectors u , h , and y , respectively.To succinctly represent the convolution, let T ∈ R m×n be the Toeplitz matrix formed by the input u , allowing us to express h ⊛ u = hT ⊤ .Further- more, assuming x ∈ R n to be an impulse response variable, we introduce X ∈ R n×n as a Hankel matrix formed by the entries of x .From [58-61], the minimum order time domain system identification problem can be formulated as: (60b) ensures that X is a Hankel matrix and (60c) holds to make the result by applying the input, u , to the optimal impulse response, x , fit the available noisy data, ŷ , in a non- trivial sense.Defining the convex set can be cast as: which is clearly identical to the problem in (1).The problem was solved using the same pQN-ADMM approach discussed in Sect. 4. Let T = m = 50 and n = 40 .It is pertinent to note that m < T + n − 1 , is a reasonable assumption aligning with practical applications where only a specific window is available to observe the output.The simulation is conducted across 5 distinct original system orders denoted by η ∈ {2, 4, 6, 8, 10} .An input vector, u , is generated, with its elements being independent and following a uniform distribution over the interval [−5, 5] .For each η :

Rank(X),
1 Fifty random stable systems are generated using the 'drss' command in MATLAB. 2 The generated input is applied to each system, yielding the corresponding noisy output ŷ. 3 Given the output ŷ , the problem specified in ( 60) is solved, and the rank of the cor- responding system is computed using singular value decomposition. 4 The obtained results are averaged to derive the corresponding average rank for each original η.
Figure 3 presents the average rank results obtained through the nuclear norm and pQN-ADMM heuristics.The outcomes correspond to two distinct threshold values, wherein the threshold is defined as the value below which the singular value is considered zero.Notably, the introduced pQN-ADMM approach demonstrates superior performance compared to the nuclear norm heuristic for both threshold values.Furthermore, as the threshold value decreases from 10 −4 to 10 −5 , the pQN-ADMM's behavior remains consistent, while the average rank for the nuclear norm exhibits an increase.This observation underscores the robustness of the derived pQN-ADMM relative to the nuclear norm approach.
Table 1 provides the standard deviation values for the algorithms.It is evident that the standard deviation remains constant for the pQN-ADMM when altering the threshold; conversely, it increases for the nuclear norm as the threshold value decreases.

Accelerated proximal gradient (APG) algorithm
In this section, we present numerical results for the APG method, as outlined in Algorithm 3. Our primary objective is to address the minimization problem (45) with f (x) = �Ax − b� 2 .Consistent with the approach in [62], we initiate the process by generating the target signal x * through: where the design parameters ⊂ [n] , and � (1)  i , � (2)  i for i ∈ are chosen as follows: 1 The index set ⊂ [n] is constructed by selecting a subset of [n] with cardinality s uniformly at random; 2 {� (1)  i } i∈� are Independent and Identically Distributed (IID) Bernoulli random vari- ables taking values ±1 with equal probability; 3 {� (2)  i } i∈� are IID uniform [0, 1] random variables.
The measurement matrix A ∈ R m×n is constructed as a partial Discrete Cosine Trans- form (DCT) matrix, with its rows corresponding to m < n frequencies.Specifically, these m indices are selected uniformly at random from the set [n].The noisy measurement vector b ∈ R m is subsequently defined as b = A(x * + ǫ 1 ) + ǫ 2 , where ǫ 1 and ǫ 2 are IID random vectors with entries following zero mean Gaussian distributions with variances σ 2 1 and σ 2 2 , respectively.In our experiments, n = 4096 , s = ⌈0.5m⌉and the APG algorithm memory to 5, i.e., l = 5 in Algorithm 3. Following the medium noise setup in [63], we set σ 1 = 0.005 , σ 2 = 0.001.
For the objective function f (x) = �Ax − b� 2 , the Lipschitz constant is given by L = 2�A� 2 .Our experimental design encompasses varying values of m, representing the number of noisy measurements, and µ , serving as the trade-off parameter in (45).For each unique combination of (m, µ) , we conduct 20 random instances of the triplet (x * , A, b) to account for the inherent statistical variability of the problem.Each ran- dom instance is subsequently solved using Algorithm 3, and the average performance is reported.The termination criterion for Algorithm 3 is defined as the relative error between consecutive iterates satisfies x k − x k−1 / x k−1 ≤ 10 −5 .
In our experiments, we conducted a comparative analysis of solving (45) for p = 0.5 against p = 1 , corresponding to ℓ 1 -optimization for sparse recovery.Specifically, for p = 0.5 , denoting ℓ 0.5 minimization, we employed Algorithm 3, referred to as ℓ 0.5 exact.Additionally, we utilized Algorithm 2 from [64], denoted as ℓ 0.5 approx.Conversely, for p = 1 , where the ℓ 1 -minimization problem is convex, we employed the FISTA algorithm from [42].The solutions are denoted as x , while the target signal, derived from (62), is denoted as x * .In Algorithm 3, we initialized x 0 as a zero vector, and x 1 was set to the ℓ 1 norm solution.
Figures 4 and 5 illustrate the relationship between average error, sparsity, and µ for var- ious values of n/m.A discernible trend is observed wherein the average error decreases while sparsity increases with an increase in µ .When µ is small, greater emphasis is (62) x * i = � (1)  placed on loss function, emphasizing ℓ 0.5 quasi-norm minimization.Consequently, the sparsity level, as depicted in Fig. 5, remains low.Conversely, for higher values of µ , more weight is assigned to the regularization term's minimization, resolving �Ax b� 2 , resulting in decreased error (Fig. 4) accompanied by an increase in sparsity.
Figure 6 provides insight into the statistics of the number of iterations until convergence for both the ℓ 0.5 exact and approximate algorithms.Notably, with a sufficient number of available realizations, specifically for n/m = 8 and n/m = 16 , both algo- rithms require approximately the same number of iterations.However, as the number of available realizations decreases, particularly for n/m = 32 and higher, the exact proximal solution exhibits a significantly lower number of iterations to converge.This observation, coupled with the findings in Figs. 4 and 5, suggests that our algorithm not only yields a comparable solution to the approximate method but also converges with fewer iterations.

Conclusion
In this paper, we introduced a non-convex ADMM algorithm, denoted as pQN-ADMM, designed for solving the ℓ p quasi-norm minimization problem.Significantly, our proposed algorithm serves as a versatile approach for tackling ℓ p problems, as it does not rely on specific structural assumptions for the convex constraint set.Moreover, we delved into the problem of solving a non-convex relaxation of RMPs utilizing the Schatten-p quasi-norm.This relaxation was established as the ℓ p minimization of the singular of the variable matrix, rendering it amenable to the pQN-ADMM algorithm.For scenarios involving constraints defined by differentiable functions with Lipschitz continuous gradients, a proximal gradient step was employed, mitigating the need for a convex projection step.This enhancement not only accelerates the algorithm but also ensures its convergence.Illustrating the numerical results, we applied the pQN-ADMM to diverse examples, encompassing sparse vector reconstruction, matrix completion, and system identification.The algorithm demonstrated competitiveness against variousℓ p -based baselines, underscoring its efficacy across a spectrum of applications.
the problem in(29) can be written as:To formulate the problem in a manner amenable to ADMM, we introduce auxiliary variables, Y ∈ R m×n and z = [z i ] i∈[L]  .This transformation leads to the following representa- tion of the problem in (30): ), we evaluate the average of 50 runs based on two metrics: a) the Relative Frobenius Distance (RFD) to the matrix M , defined as RFD = �X * −M� f �M� f , and b) the Relative Error to Singular (REtS) values of M , given by REtS i = |σ i (X * )−σ i (M)| σ i (M) for i ∈ [min{m, n}].

Fig. 5
Fig. 5 Sparsity µ for different values of n/m.Yellow and cyan shades are the standard deviations for the exact and approximate ℓ 0.5 quasi-norms, respectively

Fig. 4
Fig. 4 Average error vs µ for different values of n/m.Yellow and cyan shades are the standard deviations for the exact and approximate ℓ 0.5 quasi-norms, respectively

1
Standard deviation for different threshold values Average rank vs. original system order.Red and cyan colors are for the nuclear norm and pQN-ADMM algorithm, respectively