Accelerated Proximal Dogleg Majorization for Sparse Regularized Quadratic Optimization Problem

TL;DR

This paper introduces PDOME, an accelerated algorithm for sparse regularized quadratic optimization that improves convergence speed and guarantees global convergence to critical points, surpassing the original PDOM method.

Contribution

The paper proposes PDOME, a novel accelerated algorithm with extrapolation and enhanced line search, providing global convergence analysis for nonconvex nonsmooth problems.

Findings

PDOME converges faster than PDOM in numerical experiments.

PDOME achieves better local optima in practice.

Theoretical analysis confirms global convergence under KL property.

Abstract

This paper addresses the problems of minimizing the sum of a quadratic function and a proximal-friendly nonconvex nonsmooth function. While the existing Proximal Dogleg Opportunistic Majorization (PDOM) algorithm for these problems offers computational efficiency by minimizing opportunistic majorization subproblems along mixed Newton directions and requiring only a single Hessian inversion, its convergence rate is limited due to the nonconvex nonsmooth regularization term, and its theoretical analysis is restricted to local convergence. To overcome these limitations, we firstly propose a novel algorithm named PDOM with extrapolation (PDOME). Its core innovations lie in two key aspects: (1) the integration of an extrapolation strategy into the construction of the hybrid Newton direction, and (2) the enhancement of the line search mechanism. Furthermore, we establish the global convergence of the entire sequence generated by PDOME to a critical point and derive its convergence rate under the Kurdyka-Lojasiewicz (KL) property. Numerical experiments demonstrate that PDOME achieves faster convergence and tends to converge to a better local optimum compared to the original PDOM.

Tables

Table 1. Table 1 : Summarize whether the existing methods incorporate second-order information (SOI), adopt opportunistic majorization (OM), utilize extrapolation techniques, compute the inverse of the Hessian matrix once, and make assumption on r r , as well as the convergence of subsequences 3 3 3 “-” means not given. .

Algorithm	Convexity of $r$	SOI	Extrapolation	OM Strategy	1 Hessian Inversion	Convergence
PG [1]	Yes	No	No	No	No	Global
PN [12]	Yes	Yes	Yes	No	No	Local
mAPG [13]	Yes	No	Yes	No	No	Global
APGnc [14]	Yes	No	Yes	No	No	Local
PANOC [15]	No	Yes	No	No	No	Global
PANOC+ [16]	No	Yes	Yes	No	No	Global
PDOM [17]	No	Yes	No	Yes	Yes	Local
sPDOME (ours)	No	Yes	Yes	Yes	Yes	-
PDOME (ours)	No	Yes	Yes	Yes	Yes	Global

Table 2. Table 2 : Average NRE and number of iterations to reach ‖ u k + 1 ‖ < 10 − 12 \|u^{k+1}\|<10^{-12} , u k + 1 ∈ ∂ Q ​ ( x k + 1 ) u^{k+1}\in\partial Q(x^{k+1}) for 20 independent trials with different m m . Sparsity level = 0.01 m m .

$m$	100	500	1000
Algorithm	NRE/#Iter	NRE/#Iter	NRE/#Iter
PG	7.00e-01 / 140.8	1.82e-01 / 171.2	9.24e-02 / 117.8
mAPG	1.00e+00 / 2.0	4.88e-01 / 182.1	2.91e-01 / 134.1
PANOC	3.41e+01 / 29.2	6.24e-01 / 38.8	4.56e-01 / 5.2
PANOC+	8.12e-14 / 13.4	1.30e-02 / 36.0	2.29e-13 / 44.0
PDOM	5.10e-13 / 151.8	2.40e-13 / 38.0	1.74e-13 / 27.2
sPDOME	8.68e-13 / 15.4	6.12e-13 / 18.9	6.88e-13 / 18.7
PDOME	9.11e-13 / 16.7	7.26e-13 / 24.9	6.20e-13 / 27.3

Table 3. Table 3 : Average NRE and number of iterations to reach ‖ u k + 1 ‖ < 10 − 12 \|u^{k+1}\|<10^{-12} , u k + 1 ∈ ∂ Q ​ ( x k + 1 ) u^{k+1}\in\partial Q(x^{k+1}) for 20 independent trials with different m m . Sparsity level = 0.01 m m .

$m$	100	500	1000
Algorithm	NRE/#Iter	NRE/#Iter	NRE/#Iter
PG	2.36e-02 / 348.3	5.25e-02 / 122.7	5.24e-02 / 87.3
mAPG	3.08e-05 / 640.5	3.92e-05 / 243.9	4.11e-05 / 102.7
PANOC	1.19e+00 / 971.9	1.17e-03 / 122.8	1.26e-03 / 5.0
PANOC+	2.36e-04 / 296.5	3.59e-04 / 123.5	3.89e-04 / 56.6
PDOM	7.13e-04 / 1621.1	1.08e-03 / 55.8	1.17e-03 / 232.8
sPDOME	3.12e-04 / 426.1	1.90e-04 / 138.6	1.81e-04 / 47.5
PDOME	1.24e-04 / 1377.1	1.34e-04 / 313.7	1.36e-04 / 97.6

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }Q(x):=s(x)+r(x),$$

$$\underset{x\to x_0}{\lim }\inf Q(x)\ge Q\left(x_0\right).$$

$$\Vert \nabla f(x)-\nabla f(y)\Vert \le L_f\Vert x-y\Vert .$$

$$\underset{y\to x,y\ne x}{\mathrm{liminf}}\frac{f(y)-f(x)-⟨v,y-x⟩}{\Vert y-x\Vert }\ge 0.$$

$$\partial f(x):=\{v\in {\mathbb{R}}^n:\exists x^k\to x,f\left(x^k\right)\to f(x),v^k\in \hat{\partial }f\left(x^k\right)\to v,k\to \infty \}.$$

$${\psi }^{\prime }(r(x)-r(\hat{x}))\mathrm{dist}(0,\partial r(x))\ge 1.$$

$$\mathrm{dist}(0,\partial r(x))\ge C(r(x)-r(\hat{x}){)}^{\theta },$$

$${\varphi }^{\prime }(h(z)-h(\bar{z}))\text{dist}(0,\partial h(z))\ge 1,$$

$$\underset{x}{\min }\underset{m_{pg}(x;y^k)}{\underbrace{s\left(y^k\right)+⟨\nabla s\left(y^k\right),x-y^k⟩+\frac{1}{2\eta }{\Vert x-y^k\Vert }^2}}+r(x)$$

$$=\frac{1}{2\eta }{\Vert x-\left(y^k-\eta \nabla s\left(y^k\right)\right)\Vert }^2+r(x)+c,$$

$$y^k=x^k+\zeta \left(x^k-x^{k-1}\right).$$

$$Q(x^{k+1})$$

$$\le m_{pg}(x^{k+1};y^k)+r(x^{k+1})$$

$$\le m_{pg}(x^k;y^k)+r(x^k),$$

$$\displaystyle d(\mu):=\left\{\begin{array}[]{ll}d_{\eta}&\mu\in(0,1],\\ d_{\eta}+(\mu-1)(d_{N}-d_{\eta})&\mu\in(1,2],\end{array}\right.$$

$${\lambda }_{\min }\Vert g{\Vert }^2\le \Vert g{\Vert }_{\mathrm{M}}^2=g^T\mathrm{M}g\le {\lambda }_{\max }\Vert g{\Vert }^2,$$

$$s(x):=⟨g,x⟩+\frac{1}{2}\Vert x{\Vert }_{\mathrm{M}}^2,m_{\mu }(x;y^k):=⟨g_{\mu },x⟩+\frac{1}{2{\eta }_{\mu }}\Vert x{\Vert }^2.$$

$$⟨d(\mu ),\nabla s(d(\mu ))⟩\le 0,$$

$$⟨d(\mu ),\nabla s(d(\mu ))⟩$$

$$={\eta }^2\left(\Vert g{\Vert }_{\mathrm{M}}^2-\frac{1}{\eta }\Vert g{\Vert }^2\right)$$

$$\le {\eta }^2\left({\lambda }_{\max }\Vert g{\Vert }^2-\frac{1}{\eta }\Vert g{\Vert }^2\right)$$

$$\le 0.$$

$$⟨d(\mu ),\nabla s(d(\mu ))⟩$$

$$<(\mu -1)\left(\eta \Vert g{\Vert }^2-\Vert g{\Vert }_{{\mathrm{M}}^{-1}}^2\right)$$

$$\le 0.$$

$$m_{\mu }(x;y^k):$$

$$=s(y^k)+\frac{1}{2{\eta }_{\mu }}{\Vert x-\left(y^k+d(\mu )\right)\Vert }^2-\frac{{\eta }_{\mu }}{2}{\Vert g_{\mu }\Vert }^2,$$

$${\eta }_{{\mu }^k}$$

$$=-\frac{\Vert d({\mu }^k){\Vert }^2}{({\mu }^k-2)\eta \Vert \nabla s(y^k){\Vert }^2+(1-{\mu }^k){\Vert \nabla s(y^k)\Vert }_{{\mathrm{M}}^{-1}}^2}.$$

$$x^{k+1}$$

$$=\arg \underset{x}{\min }r(x)+\frac{1}{2{\eta }_{{\mu }_k}}\Vert x-(y^k+d({\mu }^k)){\Vert }^2.$$

$$m_{\mu }(x;y^k)\ge s(x),\forall x\in \mathrm{dom}f.$$

$${\mathcal{X}}_{\mu }:=\left\{x(\beta ):=\beta d(\mu ):\forall \beta \in \mathbb{R}\right\}.$$

$${\bar{s}}^{\prime }(0)=\beta d(\mu {)}^{\mathsf{T}}\mathrm{M}d(\mu )+g^{\mathsf{T}}d(\mu ){|}_{\beta =0}=⟨g,d(\mu )⟩,$$

$${\bar{m}}_{\mu }^{\prime }(0)=g_{\mu }^{\mathsf{T}}d(\mu )+\frac{1}{{\eta }_{\mu }}\beta \Vert d(\mu ){\Vert }^2{|}_{\beta =0}=⟨g,d(\mu )⟩.$$

$${\bar{m}}_{\mu }^{\prime }(1)={g_{\mu }^{\mathsf{T}}d(\mu )+\frac{1}{{\eta }_{\mu }}\beta \Vert d(\mu ){\Vert }^2|}_{\beta =1}=0.$$

Full text

**Keywords: ** Majorization-minimization, Nonconvex and nonsmooth optimization, Proximal Newton-like method, Extrapolation, Line search.

1 Introduction

This paper investigates a class of nonconvex composite optimization problems, specifically considering objective functions formed by the sum of a convex quadratic function and a nonconvex nonsmooth function:

[TABLE]

where ${x}$ is the decision variable, $s$ is a quadratic function with $\nabla^{2}s({x})\succcurlyeq 0$111When $\nabla^{2}{s}$ is positive semi-definite, an $\epsilon\bm{I}$ can be added into $\nabla^{2}{s}$ where $\epsilon>0$ is small. and $r:\mathbb{R}^{n}\to\mathbb{R}\cup\{\infty\}$ is a nonconvex and nonsmooth function (which may represent regularization, constraints, or complex structures). We further assume that $Q$ is lower bounded, i.e., there exists a real number $h$ such that $\forall{x}\in\mathbb{R}^{n}$, $Q({x})\geq h$. Additionally, $r$ is proximal-friendly, which means that proximal operator $\operatorname{prox}_{\eta r}(\cdot):=\arg\min_{\bm{x}\in\mathbb{R}^{n}}\{r({x})+\frac{1}{2\eta}\|{x}-\cdot\|^{2}\}$ ( with step size $\eta>0$ ) is easy to compute [1].

The structured optimization problem defined in (1) arises in diverse signal processing and machine learning applications. A quintessential example is sparse signal recovery [2], which underpins techniques like channel estimation [3], audio processing [4], and blind source separation [5]. To enforce sparsity, numerous regularization strategies are employed, encompassing both convex methods, like the $l_{1}$ norm [6] and nonconvex counterparts including the $l_{1/2}$ quasi-norm [7], the $l_{0}$ pseudo-norm [8], and the capped-$l_{1}$ penalty [9]. Beyond sparse modeling, similar composite optimization formulations also appear in other key areas of machine learning, such as low-rank matrix completion [10] and robust principal component analysis (RPCA) [11].

Various algorithms exist for solving (1), among which the proximal gradient (PG) method [1] is the most widely adopted due to its simplicity and ease of implementation. PG method iteratively combines gradient descent on the smooth component with a proximal update for the nonsmooth term. However, in nonconvex settings, PG method suffers from slow convergence, exhibiting only sublinear rates of $O(1/k)$ in the worst case (where $k$ is the iteration index) [18, 19]. To accelerate convergence, techniques incorporating momentum (e.g., the accelerated proximal gradient method for nonconvex programming, APGnc [14]) or Nesterov extrapolation (e.g., the accelerated proximal gradient (APG) method [13, 20]) have been developed. These methods employ adaptive mechanisms to dynamically select between standard PG updates and accelerated variants based on objective function values. Although these accelerated variants offer improvements, their guarantees of faster convergence often require specific problem structures or assumptions [13, 18, 19]. Consequently, the desired acceleration may fail to materialize in more general, unconstrained, or challenging nonconvex optimization scenarios.

Newton-type algorithms, exemplified by the proximal Newton method (PN) [12], have garnered significant recent attention for minimizing objectives comprising a twice-differentiable term and a proximal-friendly term. At its core, each iteration involves constructing a scaled proximal operator (SPO) derived from the Hessian of the differentiable term and solving the resulting subproblem. Crucially, when the objective is convex and the SPO can be efficiently solved, this approach achieves a superlinear asymptotic convergence rate. However, the proximal Newton method faces two fundamental limitations:(1) solving the SPO presents a significant computational challenge, and the development of efficient solvers has yet to be achieved; (2) its established fast convergence is confined to convex problems, providing no assurance for the nonconvex case.

Instead of directly addressing the computationally challenging SPO, quasi-Newton approaches [15, 16, 21, 22, 23] strategically minimize the forward-backward envelope (FBE), which shares the same local minimizers as the original objective function. These methods iteratively compute the FBE’s gradient and update a quasi-Newton direction (using BFGS or L-BFGS [24]) based on this gradient, subsequently performing a line search along the minimization, it introduces two significant drawbacks: (1) the iterative updating of the Hessian approximation via BFGS/L-BFGS imposes a notable computational burden and memory overhead; (2) crucially, these approximations fail to exploit potential special structure present in the exact Hessian, potentially discarding valuable problem-specific information that could accelerate convergence.

Several algorithms leverage the core Newton-type principle of using approximate second-order information for fast convergence. Sepcially, Proximal Averaged Newton-type method for Optimal Control (PANOC) [15] innovatively integrates the forward-backward (FB) method and the FBE method to overcome their individual drawbacks, demonstrating effectiveness in nonlinear constrained optimal control and achieving superlinear convergence under mild assumptions. Building on PANOC, PANOC+ [16] addresses specific shortcomings (including those of its derivative, the PG algorithm) by introducing an adaptive step size rule tailored to the PG oracle. Validated through case studies, PANOC+ offers a complete convergence theory (handling local Lipschitz continuity) and robustness against suboptimal PG subproblem solutions. However, both PANOC and PANOC+, relying on an FBE penalty method within the Augmented Lagrangian Method (ALM) framework, can struggle with nonconvex nonsmooth constraints, often requiring approximations that degrade convergence speed. Complementing these, the PDOM algorithm [17] employs a majorization-minimization (MM) approach for problem (1), constructing a dogleg surrogate model that strategically combines gradient and Newton directions with proximal mapping. A key advantage of PDOM is that the quadratic nature of the smooth term keeps the Hessian constant, allowing its inverse to be precomputed and stored, thus eliminating the need for costly iterative matrix inversions typical of quasi-Newton methods. While PDOM provides theoretical guarantees for convergence to critical points and analyzes its local convergence rate, a significant limitation is that its convergence analysis, particularly the rate, is confined to the local domain,lacking established global convergence guarantees.

To address the limitations of the PDOM algorithm in [17], we propose a simple PDOM algorithm called sPDOME, which incorporates an extrapolation parameter mechanism to significantly improve numerical performance. To further ensure global convergence, we develop the PDOME algorithm, which innovatively combines extrapolation techniques with an improved backtracking line search, thereby establishing a rigorous theoretical framework for convergence analysis. The main contributions include:

• •

The core idea of PDOME algorithm is based on the majorization-minimization (MM) framework and incorporates extrapolation acceleration techniques. This algorithm constructs a surrogate function along the dogleg path at the extrapolated point, integrating the gradient direction and Newton-type search direction: the gradient direction ensures the reliability of the sequence of iterates, while the Newton direction accelerates the local convergence rate. Further, optimizing the line search criterion helps determine a more optimal step size, thereby supporting subsequent convergence analysis.

• •

Through theoretical analysis, we show that the limit points of the sequences generated by PDOME are critical points of the objective function. Then, by exploiting different cases of the Kurdyka-Lojasiewicz (KL) property of the objective function, we establish comprehensive convergence rate guarantees for PDOME, systematically characterizing its behavior across three distinct convergence regimes determined by the Lojasiewicz exponent. The theoretical analysis demonstrates that PDOME maintains computational efficiency comparable to conventional methods while requiring fewer proximal operator computations per iteration.

• •

We conducted numerical experiments on several well-known nonconvex and nonsmooth problems. The results show that, compared with other benchmark algorithms, both the sPDOME algorithm and the PDOME algorithm can converge quickly.

The rest of this paper is structured as follows: Section 2 presents mathematical preliminaries and related preparatory work; Section 3 introduces the proposed sPDOME and PDOME algorithms, including a detailed analysis of PDOME algorithm convergence properties; Section 4 reports relevant experimental results; Section 5 summarizes the paper and outlines future research directions.

2 Preliminaries

In this paper, $\mathbb{R}^{n}$ is defined as the $n$ dimensional Euclidean space. The symbols $\cdot,\langle\cdot,\cdot\rangle$ and $T$ represent the standard product, inner product and transpose in the space $\mathbb{R}^{n}$. For an arbitrary vector ${x}\in\mathbb{R}^{n}$, the $\ell_{2}$ norm, the $\ell_{1}$ norm, and the $\ell_{0}$ pseudo-norm are defined as $\|{x}\|:=\sqrt{{x}^{T}{x}}$, $\|x\|_{1}:=\sum_{i=1}^{n}|x_{i}|$, and $\|x\|_{0}:=|\operatorname{supp}(x)|$ where $\operatorname{supp}(\cdot)$ counts the number of nonzero elements in $x$. Given a positive semidefinite matrix $\bm{\mathrm{M}}\in\mathbb{R}^{n\times n}$, the scaled norm of $x$ is defined as $\|x\|_{\bm{\mathrm{M}}}:=\sqrt{x^{T}\bm{\mathrm{M}}x}$. Given a closed set $\Omega\subseteq\mathbb{R}^{n}$, $\operatorname{dist}(x,\Omega):=\inf\{\|y-x\|_{2}:y\in\Omega\}$ calculates the distance between $x$ and $\Omega$.

Definition 2.1**.**

(Lower semicontinuous [25]). A function $Q:\mathbb{R}^{n}\to(-\infty,+\infty]$ is said to be proper if $\operatorname{dom}Q\neq\emptyset$, where $\operatorname{dom}Q=\{x\in\mathbb{R}^{n}:Q({x})<+\infty\}$, and lower semicontinuous at point $x_{0}$ if

[TABLE]

Definition 2.2**.**

(Gradient Lipschitz continuity [25]). Let $L_{f}\geq 0$. A differentiable function $f:\mathbb{R}^{n}\to\mathbb{R}$ is said to have a Lipschitz continuity of the gradient if for all $x,{y}\in\operatorname{dom}f$ it holds that

[TABLE]

The corresponding Lipschitz constant is denoted as $L_{f}$.

The value of $L_{f}$ for a twice differentiable function $f:$ $\mathbb{R}^{n}\to\mathbb{R}$ with a positive semi-definite Hessian matrix $\bm{\mathrm{M}}\in\mathbb{R}^{n\times n}$ is the largest eigenvalue of $\bm{\mathrm{M}}$, denoted by $\lambda_{\mathrm{max}}(\bm{\mathrm{M}}).$

Definition 2.3**.**

(Subdifferential [26]). Let $f:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ be a proper and lower semicontinuous function. For a given $x\in\operatorname{dom}f$, the $Fr\acute{e}chet$ subdifferential of $f$ at $x$, written as $\hat{\partial}f({x})$, is the set of all vectors ${v}\in\mathbb{R}^{n}$ which satisfy

[TABLE]

For $x\notin\operatorname{dom}f$, we set $\hat{\partial}f({x}):=\emptyset$. The subdifferential (which is also called the limiting subdifferential) of $f$ at $x\in\mathbb{R}^{n}$, written as $\partial f(x)$, is defined by

[TABLE]

As before $\partial f(x):=\emptyset$ for $x\notin\mathrm{dom}f$, and its domain is $\mathrm{dom}~\partial f:=\{x\in\mathbb{R}^{n}:\partial f(x)\neq\emptyset\}$.

Definition 2.4**.**

*(KL property [27]). A proper closed function $r:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ is said to have the KL property at $\hat{x}\in\operatorname{dom}\partial r$ if there exists $\eta\in(0,+\infty]$, a neighborhood $\mathcal{B}_{\rho}(\hat{x})\triangleq\{x:\|x-$ $\hat{x}\|<\rho\}$, and a continuous desingularizing concave function $\psi:[0,\eta)\to[0,+\infty)$ with $\psi(0)=0$ such that

(i) $\psi$ is a continuously differentiable function with $\psi^{\prime}(x)>$ 0, $\forall x\in(0,\eta)$,

(ii) for all $x\in\mathcal{B}_{\rho}(\hat{{x}})\cap\{{u}\in\mathbb{R}^{n}:r(\hat{{x}})<r({x})<r(\hat{{x}})+\eta\},$ it holds that*

[TABLE]

A proper closed function $r$ satisfying the KL property at all points in $\operatorname{dom}\partial r$ is called a KL function.

Definition 2.5**.**

(Eojasiewicz exponent [28]). For a proper closed function $r$ satisfying the KL property at $\hat{x}\in\operatorname{dom}\hat{\partial}r$, if the desingularizing function $\psi$ can be chosen as $\psi(t)=$ $\frac{C}{1-2\theta}t^{2\theta-1}$ for some $C>0$ and $\theta\in[0,1)$, i.e., there exist $\rho>0$ and $\eta\in(0,+\infty]$ such that

[TABLE]

where $x\in\mathcal{B}_{\rho}(\hat{x})$ and $r(\hat{x})<r({x})<r(\hat{{x}})+\eta$, then we say that $r$ has the KL property at $\hat{x}$ with an exponent of $\theta.$ We say that $r$ is a KL function with an exponent of $\theta$ if $r$ has the same exponent $\theta$ at any $\hat{x}\in\operatorname{dom}\partial r,$ where the desingularizing function $\psi$ can be chosen specifically as $\psi(t)=$ $\frac{C}{\theta}t^{\theta}$ with constants $C>0$ and $\theta\in(0,1]$.

The KL property holds for a large family of functions used in optimization. For instance, all proper and closed semi-algebraic or subanalytic functions satisfy the KL property with an associated Lojasiewicz exponent $\theta\in[0,1)$ [27]. The global convergence rate of PDOME is determined by the value of the Lojasiewicz exponent. In Subsection 3.2, we provide the exponent value of the problem under test.

Lemma 2.1**.**

(Uniformized KL Property [29]). Let $\mathbb{W}$ be a compact set and $r:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{\infty\}$ be a proper and lower semicontinuous function. We assume that $r$ is constant on $\mathbb{W}$ and satisfies the KL property at each point of $\mathbb{W}$. Then, there exist $\epsilon>0$, $\eta>0$ and $\phi\in Y_{\eta}$ such that for all $\bar{z}\in\mathbb{W}$, one has

[TABLE]

for all $z\in\{z\in\mathbb{R}^{n}\mid\text{dist}(z,\mathbb{W})<\epsilon\}\cap\{z\in\mathbb{R}^{n}\mid h(\bar{z})<h(z)<h(\bar{z})+\rho\}$.

2.1 PG Method from the Majorization-minimization Angle

In this subsection, we examine the proximal gradient method from the perspective of a majorization-minimization algorithm and identify certain limitations that lead to its slow convergence rate. The PG method is a well-established method for solving composite optimization problems of the form (1). At each iteration, the PG method performs a proximal line search in the direction of the negative gradient, employing a positive step size $\eta$. At a given point $y^{k}$, with $\zeta$ as the momentum coefficient and $c$ being a constant, PG method solves the following surrogate function

[TABLE]

where

[TABLE]

Assuming that $r(\cdot)$ is proximable, the solution of (7) can be obtained efficiently from a computational standpoint. The iterates generated by the algorithm yield a nonincreasing sequence of objective function’s values. This property is ensured by the fact that the surrogate function $m_{pg}({x};{y}^{k})$ serves as an upper bound for $s({x})$. Specifically, $m_{pg}({x};{y}^{k})\geq s({x})$ holds for all ${x}\in\operatorname{dom}Q$, provided that $\eta<1/L_{s}$.

[TABLE]

where the first inequality follows from the majorization step, and the second arises from the proximal operator’s property. A proximal gradient method with extrapolation and line search (PGels) is proposed in [30] to address composite optimization problems that are potentially nonconvex, nonsmooth, and non-Lipschitz. By constructing an auxiliary function, the global subsequential convergence of PGels is proved. With appropriate parameter selection, PGels can be simplified to PG and PGe. The convergence guarantee of $Q\left(x^{k}\right)$ can be treated in the same way as in [30].

However, using the negative gradient direction often results in slow convergence, especially for nonconvex functions [9, 19, 31]. Furthermore, in cases where the Hessian matrix exhibits a large condition number, gradient-based methods become inefficient because of excessively slow convergence, as noted in [32, Chapter 9].

2.2 Hybrid Direction and Opportunistic Majorization

To address the limitations identified in Subsection 2.1, we propose the sPDOME and PDOME algorithm. These approaches integrates a dogleg search strategy, drawing inspiration from trust region methods. At each iteration, the algorithm constructs and minimizes an opportunistically majorized surrogate function along the dogleg path, replacing conventional gradient-based updates.

Given $\mu\in(0,2]$, the dogleg path is denoted as

[TABLE]

where $g:=\nabla s(y)$, $d_{\eta}:=-\eta\nabla s(y)$, and $d_{N}:=-(\nabla^{2}s(y))^{-1}\nabla s(y)$, $\eta$ is the fixed step size of the gradient direction, with $\eta\in(0,1/L_{s})$, and $d_{N}$ denotes Newton point. The gradient direction is essential for ensuring convergence to a critical point, because at $x^{\mathrm{cri}}$, the first-order optimality condition of (1) implies $0\in\partial Q(x^{\text{cri}})=\nabla s(x^{\text{cri}})+\partial r(x^{\text{cri}})$, for which the gradient direction is necessary. The construction of our path diverges fundamentally from the approach in [33, Chapter 4], specifically by excluding the scaling factor $\mu$ in the initial segment. The adoption of this modification is warranted as $d_{\eta}$ intrinsically functions as a descent direction, all while maintaining freedom from trust-region restrictions. However, the path continues to ensure descent with respect to the quadratic term.

Given a positive definite matrix $\bm{\mathrm{M}}\succ 0$ with bounded eigenvalues, it holds that

[TABLE]

where $\lambda_{\mathrm{max}}$ and $\lambda_{\mathrm{min}}$ denote the largest and the smallest eigenvalue of $\bm{\mathrm{M}}$, respectively. For the sake of subsequent analysis, we perform a coordinate transformation that shifts the origin to the point $(y^{k},s(y^{k}))$. Under this new coordinate system, the smooth part of the objective function and its corresponding surrogate are rewritten as

[TABLE]

3 Main Results

In this section, we introduce the PDOME algorithm for solving problem (1) which may involve nonconvexity and nonsmoothness, aiming to address the limitations outlined in Subsection 2.1. Building on the existing technical framework, we further propose a novel extrapolation technique to enhance the algorithm’s performance. Specifically, on the basis of extrapolation, we integrate the dogleg search strategy originally developed in the trust region methodology into the PDOME algorithm, replacing the conventional gradient-based descent direction with a hybrid search direction. This integration enables the algorithm to utilize more sophisticated directional information, thereby potentially accelerating convergence and improving robustness.

Lemma 3.1**.**

The inequality in the following equation,

[TABLE]

is satisfied when $\eta=\frac{1}{\lambda_{\max}}$, where $\lambda_{\mathrm{max}}$ denote the largest eigenvalue of $\bm{\mathrm{M}}$. For any other $\eta\in(0,\frac{1}{\lambda_{\max}})$, the strict inequality holds.

Proof.

We begin with the trivial case where $\mu\in(0,1]$,

[TABLE]

Then for the remaining range of $\mu\in(1,2]$, we have

[TABLE]

The equality of (16) only holds when $\eta$ is exactly $1/\lambda_{\max}$. ∎

The effectiveness of an MM algorithm fundamentally depends on the construction of a surrogate function that acts as a tight upper bound for the original objective function. In PDOME, the local surrogate function $m_{\mu}$ is the projection of $m_{pg}$ onto the path direction, that is

[TABLE]

where $g_{\mu}=\frac{\langle\nabla s\left(y^{k}\right),d(\mu)\rangle}{\|d(\mu)\|^{2}}d(\mu)$ and $\eta_{\mu}=-\frac{\|d(\mu)\|^{2}}{\langle\nabla s(y^{k}),d(\mu)\rangle}$, for $\mu\in(0,2]$. The step size $\eta_{\mu}$ is allowed to surpass $\eta$.

Lemma 3.2**.**

$\eta_{\mu}$* is an increasing function of $\mu\in[0,2]$.*

Proof.

The step size $\eta_{\mu}$ is considered an increasing function of $\mu$ within the interval $[0,2]$ provided that $\frac{d}{d\mu}\eta_{\mu}\geq 0$. The positive gradient can be proved with simple algebra. ∎

Lemma 3.3**.**

The sequence $\left\{\eta_{\mu^{k}}\right\}_{k\in\mathbb{N}}$ is bounded.

Proof.

Based on the definition of $\eta_{\mu}$, it holds that

[TABLE]

Since the eigenvalues of ${\bm{\mathrm{M}}}$ are bounded, both the numerator and denominator are constrained by finite values. Consequently, the sequence $\{\eta_{\mu^{k}}\}_{k\in\mathbb{N}}$ is established as bounded. ∎

In each iteration, the update rule is

[TABLE]

It is not guaranteed that the new iterate $x^{k+1}$ results in a lower objective function value, since there is no assurance that $m_{\mu}(x;y^{k})$ majorizes $s(x)$ for all choices of $\mu$. In what follows, we examine the MM condition under different ranges of $\mu$. We begin with the case where $\mu\in(0,1]$, in which $m_{\mu}$ reduces to $m_{pg}$ with $g_{\mu}$ becoming $g$ and $\eta_{\mu}$ simplifying to $\eta$. According to (7), the surrogate function $m_{\mu}$ provides a uniform upper bound on $s$, meaning that

[TABLE]

Next, we consider the case where $\mu\in(1,2]$. In this range, the surrogate function $m_{\mu}$ only majorizes $s$ along the specific line segment that connects the current iterate and the path point. Unlike the classical MM principle cited in [34, 35], where the surrogate is required to upper bound the objective function globally or over the entire domain, $m_{\mu}$ does not necessarily remain above function $s$ everywhere. We refer to this more flexible condition as OM.

Theorem 3.1**.**

For any given $\mu\in(1,2]$, consider the line connecting 0 and $d(\mu)$ which is given by

[TABLE]

Define $\bar{s}(\beta):=s({x}(\beta))$ and $\bar{m}_{\mu}(\beta):=m_{\mu}({x}(\beta);y^{k})$. It holds that $\bar{s}(\beta)\leq\bar{m}_{\mu}(\beta)$ for all $\beta\in\mathbb{R}$, or equivalently, $s(x)\leq m_{\mu}({x};y^{k})$ for all $x\in\mathcal{X}_{\mu}$.

Proof.

The proof of the theorem can be established by considering Lemma 3.4, Lemma 3.5, and Lemma 3.6. ∎

Lemma 3.4**.**

Given $\bar{s}$ and $\bar{m}_{\mu}$ defined in Theorem 3.1, it holds that $\bar{s}(0)=\bar{m}_{\mu}(0)=0$ and ${\bar{s}^{\prime}\left(0\right)}={\bar{m}_{\mu}}^{\prime}\left(0\right)<0$.

Proof.

It is easy to show that $\bar{s}(0)=s(x(0))=0$, and $\bar{m}_{\mu}(0)=m_{\mu}(x(0);y^{k})=0$. Then we prove the negative gradient. It holds that

[TABLE]

[TABLE]

From Lemma 3.1, we prove that $\bar{s}^{\prime}(0)=\bar{m}_{\mu}^{\prime}(0)<0.$ ∎

Lemma 3.5**.**

*Let $m(x;y^{k})$ and $s(x)$ be univariate strictly convex quadratic functions. Suppose that

(i) $m(0)=s(0)$ and $m^{\prime}(0)=s^{\prime}(0)\neq 0$;

(ii) $x_{m}^{\star}=\eta x_{s}^{\star}$ for some $\eta\in(0,1)$, where $x_{m}^{\star}:=\operatorname*{arg\,min}_{x}m(x)~and~x_{s}^{\star}:=\operatorname*{arg\,min}_{x}s(x),$*

then, it holds that

• •

$\begin{aligned} \left|\frac{x_{m}^{\#}}{m^{\prime}(0)}\right|<\left|\frac{x_{s}^{\#}}{s^{\prime}(0)}\right|;\end{aligned}$**

• •

$\begin{aligned} m(x;y^{k})\geq s(x)\text{ for all }x,\text{ and the equality holds if and only if }x=0.\end{aligned}$**

Proof.

As both $m(x;y)$ and $s(x)$ are quadratic, one can write them as $s(x)=s(0)+s^{\prime}(0)x+\frac{1}{2\eta_{s}}x^{2}$ and $m(x;y^{k})=s(0)+s^{\prime}(0)x+\frac{1}{2\eta_{m}}x^{2}$. It is clear that $x_{q}^{\#}=-\eta_{s}s^{\prime}(0),x_{m}^{\#}=-\eta_{m}s^{\prime}(0)$. From the assumption that $x_{m}^{\#}=\arg\min_{x}m(x;y^{k})=\eta\arg\min_{x}s(x)=\eta x_{s}^{\#}$, it holds that $\left|\frac{x_{m}^{\#}}{m^{\prime}(0)}\right|=\eta_{m}=\eta\eta_{f}<\eta_{f}=\left|\frac{x_{s}^{\#}}{s^{\prime}(0)}\right|$, or equivalently, $\frac{1}{\eta_{m}}>\frac{1}{\eta_{f}}$. Therefore, $m(x;y^{k})\geq s(x)$ where the equality holds if and only if $x=0$. Both claims in the lemma are therefore proved. Based on Lemma 3.4 and Lemma 3.5, Theorem 3.1 can be proved by showing the lemma below. ∎

Lemma 3.6**.**

Considering $\bar{s}(\beta)$ and $\bar{m}(\beta)$ defined in Theorem 3.1, it holds that $1=\arg\min_{\beta}\bar{m}_{\alpha}(\beta)\leq\arg\min_{\beta}\bar{s}(\beta).$

Proof.

It is established that

[TABLE]

The claim that $1=\arg\min_{\beta}\bar{m}_{\mu}(\beta)$ is therefore proved. We now show that $\bar{s}^{\prime}(1)\leq 0$. It is clear that

[TABLE]

where the last inequality comes from Lemma 3.1. Combining this with Lemma 3.4 that $\bar{Q}^{\prime}(0)=\bar{m}_{\mu}^{\prime}(0)<0$, it can be concluded that $1=\arg\min_{\beta}\bar{m}_{\mu}(\beta)\leq\arg\min_{\beta}\bar{s}(\beta).$ This completes the proof. ∎

Proximal line search-based algorithms implicitly employ the principle of OM, although this principle is typically not explicitly recognized or articulated [36]. Herein, we explicitly define the OM concept and integrate it into a Newton-type optimization framework.

3.1 Algorithm Development

Theorem 3.1 indicates that, by using the non-trivial surrogate function (17), the majorization condition is satisfied as long as the new iterate lies along the specified line. This ensures that the sequence $\left\{Q\left(x^{k}\right)\right\}_{k\in\mathbb{N}}$ is monotonically decreasing. To accelerate the convergence rate of the proximal gradient method, Ochs et al. [37] introduced an inertial mechanism commonly referred to as extrapolation [30, 38] into the proximal update framework. Based on the improved algorithm presented in [17], we thus develop the sPDOME algorithm, designed for solving problem (1). The corresponding iterative scheme is formulated as follows

[TABLE]

where $d_{\gamma}(\mu)=\gamma{d}(\mu)$, and $\gamma\in(0,1)$ is a constant and typically set close to 1 in numerical experiments. The new iterate is

[TABLE]

which remains easy to solve given the assumption that the standard proximal operator is computationally simple. The overall algorithm is summarized below. To find the largest $\mu^{k}$, Line 5 uses a strategy similar to that presented in [15].

Remark 3.1**.**

If $\zeta=0$, the sPDOME Algorithm 1 reduces to $PDOM$ Algorithm [17].

Building upon Algorithm 1, we made some minor adjustments to the range of the inertial coefficient $\zeta$ and improved the line search step for $\mu^{k}$ to be found, such as $\langle-\nabla s\left({y}^{k}\right)+\frac{\langle\nabla s\left({y}^{k}\right),{d}(\mu^{k})\rangle}{\|{d}(\mu^{k})\|^{2}}{d}(\mu^{k}),x^{k}-y^{k}\rangle\leq 0$, thus obtaining the PDOME method, see Algorithm 2 for details.

The PDOME algorithm is designed to terminate upon approaching a critical point $x^{\ast}$, at which the condition $0\in\partial Q(x^{\star})$ is satisfied. Based on the optimality condition associated with the proxima operator in (21), the following relation holds:

[TABLE]

this implies

[TABLE]

PDOME terminates when $\|\partial Q(\bm{x}^{k+1})\|$ is sufficiently small:

[TABLE]

where $n$ is the dimension of $x$, $\epsilon^{abs}>0$ and $\epsilon^{rel}>0$ are two small positive constants (motivated by  [39, Section 3.3]).

This stopping criterion is different from directly using $\|x^{k+1}-x^{k}\|$, commonly adopted for proximal algorithms [12]. The relationship between these two different stopping criteria can be roughly quantified by the triangle inequality

[TABLE]

As $1/\eta_{\mu^{k}}$ in (25) can be very large, small value of $\|x^{k+1}-x^{k}\|$ does not necessarily imply getting close to a critical point. Now we formally present the PDOME in Algorithm 2. To track the largest $\mu^{k}$, we employ a similarly straightforward strategy as presented in [15].

3.2 The Convergence and Convergence Rate Analysis

In this subsection, we analyze the convergence behavior of the sequence generated by the Algorithm 2, showing that it converges to a critical point of $H_{\delta}(x_{k})$. Additionally, under the KL condition, we derive the global convergence rate. We begin by establishing the monotonic decrease of the objective function values throughout the iterations.

Theorem 3.2**.**

The sequence $\left\{Q(x^{k})\right\}_{k\in\mathbb{N}}$ generated by Algorithm 2 satisfies the following inequality.

Proof.

It holds that

[TABLE]

where the first inequality follows from the backtracking rule, the second inequality holds by virtue of $0<\gamma<1$, and the third inequality from the proximal operator. ∎

Lemma 3.7**.**

*Suppose that $\left\{x^{k}\right\}_{k\in\mathbb{N}}$ is a sequence generated by Algorithm 2, then it holds that

(i) The sequence $\{H_{\delta_{k}}(x^{k})\}$ is monotonically nonincreasing. In particular, for any $k\in\mathbb{N}$, it holds*

that

[TABLE]

(ii)  $\lim_{k\to\infty}\left\|{x}^{k+1}-{x}^{k}\right\|^{2}\to 0.$

Proof.

(i) To simplify the notations in our analysis, we denote

[TABLE]

In the following, we show that the sequence is monotonically nonincreasing. By following (21), the path search procedure finds a new update ${x}^{k+1}$ to make $m_{\mu^{k}}({x}^{k+1};{y}^{k})$ an upper bound of $s({x}^{k+1})$, thus we have

[TABLE]

From equation (28), we obtain

[TABLE]

[TABLE]

Based on the convexity of the quadratic function $s$, we obtain that

[TABLE]

Using equations (30) and (31), we derive

[TABLE]

Combining (29) and (32), we derive

[TABLE]

Due to $\left\langle{-g+{g}_{\mu^{k}}},{x}^{k}-{y}^{k}\right\rangle\leq 0$ in Algorithm 2, where $g=\nabla s\left({y}^{k}\right)$ and ${g}_{\mu^{k}}=\frac{\langle\nabla s\left({y}^{k}\right),{d}(\mu^{k})\rangle}{\|{d}(\mu^{k})\|^{2}}{d}(\mu^{k})$, its corresponding geometric interpretation is illustrated in Figure 1.

[TABLE]

we know that

[TABLE]

where the inequality follows from the fact that

[TABLE]

In the second inequality of equation (34) is obtained by combining equations (8), (35), and (36). Based on the definition in equation (27) and the result in equation (34), we can derive equation (37).

[TABLE]

Due to $0<\zeta<\frac{1-\gamma}{2-\gamma}$ in Algorithm 2, then the sequence follows $H_{\delta_{k}}(x^{k})$ is monotonically nonincreasing.

(ii) Then, summing up (37) from $k=0,1,\ldots,N$ and $x^{-1}=x^{0}$, it yields

[TABLE]

Given that $\{\eta_{\mu^{k}}\}_{k\in\mathbb{N}}$ is bounded, we derive that

[TABLE]

Since $x^{k+1}$ is set to $v^{k+1}$ whenever $v^{k+1}$ yields a smaller objective value, the following condition must also be satisfied:

[TABLE]

For the proof of this statement, refer to [40]. This concludes the proof. ∎

Theorem 3.3**.**

Suppose that $H_{\delta}$ is lower-bounded, $s$ is a quadratic function, and $r$ is a lower semicontinuous function. Let $\{x^{k}\}_{k\in N}$ be a sequence generated by the Algorithm 2 converging to $x^{\star}$. Then $0\in\partial H_{\delta}(x^{\star})$, i.e., $x^{\star}$ is a critical point.

Proof.

The theorem can be proved by combining Lemma 3.3 (which establishes step-size boundedness) with part (ii) of Lemma 3.7, following the methodology of [13, Theorem 1]. ∎

We now analyze the global convergence rate of the PDOME using the KL property. Initially, we show that $g_{\mu^{k}}$ converges to the gradient direction of $s$ as $k\rightarrow\infty$.

Lemma 3.8**.**

Let $e^{k}={g}_{\mu^{k}}-\nabla s({y}^{k})$. The sequence $\left\{\left\|e^{k}\right\|\right\}_{k\in\mathbb{N}}$ converges to 0 as $k\to\infty$ and we have

[TABLE]

Proof.

By analyzing the optimality condition of (21), we obtain that

[TABLE]

By triangle inequality and smoothness of $\nabla s$, we have

[TABLE]

Then, we have

[TABLE]

Based on $\lim_{k\to\infty}\left\|x^{k+1}-x^{k}\right\|^{2}\to 0$ and $\lim_{k\to\infty}\left\|x^{k}-x^{k-1}\right\|^{2}\to 0$ in Theorem 3.3, we have $\lim_{k\to\infty}\left\|e^{k}\right\|\to 0.$ ∎

Theorem 3.4**.**

Suppose that $H_{\delta}$ satisfies the KL property on $\omega(x^{k})$ which is the cluster point set of $\{x^{k}\}_{k\in\mathbb{N}}$, then the sequence $\{x^{k}\}_{k\in\mathbb{N}}$ generated by Algorithm 2 has summable residuals, $\sum_{k=0}^{\infty}\|x^{k+1}-x^{k}\|<\infty$.

Proof.

Following the same procedure as in [19, Theorem 2.9], and considering the descent property in (37) together with the property from (41) that $g_{\mu^{k}}$ converges to the gradient direction of $s$ as $k\rightarrow\infty$. one can easily show that the sequence $\left\{{x}^{k}\right\}_{k\in\mathbb{N}}$ has a finite length. ∎

Whenever the KL property is invoked, we shall adopt the results given in Theorem 3.4. We consider a sequence $\{x^{k}\}_{k\in\mathbb{N}}$ in Algorithm 2, computed by means of an abstract algorithm satisfying the following inequality:

$\mathbf{H}_{1}$ (Sufficient decrease): From equation (37), for each $k\in\mathbb{N}$, where $\frac{\zeta(\gamma-2)+1-\gamma}{2\gamma\eta_{\mu^{k}}}>0,$

[TABLE]

$\mathbf{H}_{2}$ (Relative error): For each $k\in\mathbb{N}$, where $\frac{\gamma\eta_{u^{k}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k}}+1)+\zeta}>0$ and $\varepsilon_{k+1}=\frac{\gamma\eta_{u^{k}}\|e^{k}\|}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k}}+1)+\zeta}\geq 0$,

[TABLE]

Let us assume that $\|x^{k}-x^{k-1}\|\leq\Re\|x^{k+1}-x^{k}\|$ holds for some $\Re>0$. In conjunction with formula (42), we arrive at $\mathbf{H}_{2}$.

$\mathbf{H}_{3}$ (Parameters): The sequences $(\frac{\zeta(\gamma-2)+1-\gamma}{2\gamma\eta_{\mu^{k}}})_{k\in\mathbb{N}}$, $(\frac{\gamma\eta_{u^{k-1}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k-1}}+1)+\zeta})_{k\in\mathbb{N}}$ and $(\varepsilon_{k})_{k\in\mathbb{N}}$ satisfy

(i) $\frac{\zeta(\gamma-2)+1-\gamma}{2\gamma\eta_{\mu^{k}}}\geq\underline{a}>0$ for all $k\geq 0$;

(ii) $(\frac{\gamma\eta_{u^{k-1}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k-1}}+1)+\zeta})_{k\in\mathbb{N}}\notin l^{1}$;

(iii)$\sup_{k\in\mathbb{N}^{\star}}\frac{2\eta_{{u}^{k}}[(\Re\zeta+1)(Ls\cdot\gamma\cdot\eta_{{u}^{k-1}}+1)+\zeta]}{[\zeta(\gamma-2)+(1-\gamma)]\eta_{{u}^{k-1}}}<+\infty$;

(iv) $(\varepsilon_{k})_{k\in\mathbb{N}}\in l^{1}$.

Theorem 3.5**.**

Let $H_{\delta}$ have the KL property at a global minimum $x^{\star}$ of $H$. Let $\left\{{x}^{k}\right\}_{k\in\mathbb{N}}$ be a sequence satisfying $\mathbf{H}_{1},\mathbf{H}_{2}$ and $\mathbf{H}_{3}$ with $\epsilon_{k}\equiv 0.$ There, exist $\rho>0$ such that if $x^{0}\in\mathcal{B}_{\rho}(\hat{x})$, then $\left\{{x}^{k}\right\}_{k\in\mathbb{N}}$ has finite length and converges to a global minimum $x^{\star}$ of $H_{\delta}$.

Proof.

As noted in [19], Theorem 3.5 allows for a more general formulation. For example, when $x^{\star}$ is a local minimum of $H_{\delta}$, a growth property holds locally (refer to [19, Remark 2.11] for details).

The proof of Theorem 3.1 and Theorem 3.2 draws on the reasoning presented in [19, Section 2.3], with adaptations made to account for the existence of errors and the variability shown by the parameters. ∎

Theorem 3.6**.**

*Let $\{x^{k}\}_{k\in\mathbb{N}}$ be any sequence generated by Algorithm 2. Suppose that $H_{\delta}$ satisfies the KL property on the cluster points of $\{x^{k}\}_{k\in\mathbb{N}}$ with exponent $\theta\in(0,1)$, then $\{x^{k}\}_{k\in\mathbb{N}}$ converges to $x^{\star}$ such that $0\in\partial H_{\delta}(x^{\star})$ and the following inequalities hold. Assume $\varphi(t)=\frac{C}{\theta}t^{\theta}$ for some $C>0$, $\theta\in[0,1]$.

(i)If $\theta=1$ and $\inf_{k\in\mathbb{N}}\frac{[\zeta(\gamma-2)+(1-\gamma)]\gamma\eta_{u^{k}}}{2[(\Re\zeta+1)(Ls\cdot\gamma\cdot\eta_{u^{k}}+1)+\zeta]^{2}}>0$, then $x^{k}$ converges in finite time.

(ii)If $\theta\in[\frac{1}{2},1]$, $\operatorname*{sup}_{k\in\mathbb{N}}\frac{\gamma\eta_{u^{k-1}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k-1}}+1)+\zeta}<+\infty$ and $\operatorname*{inf}_{k\in\mathbb{N}}\frac{\zeta(\gamma-2)+1-\gamma}{2[(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k}}+1)+\zeta]}>0$, there exist $c>0$ and $k_{0}\in\mathbb{N}$ such that*

• •

$\begin{aligned} H_{\delta}(x^{k})-H_{\delta}(x^{\star})&=O\left(\exp\left(-c\sum_{n=k_{0}}^{k-1}\frac{\gamma\eta_{u^{n}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{n}}+1)+\zeta}\right)\right),\end{aligned}$**

• •

$\begin{aligned} \|x^{\star}-x^{k}\|&=O\left(\exp\left(-\frac{c}{2}\sum_{n=k_{0}}^{k-2}\frac{\gamma\eta_{u^{n}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{n}}+1)+\zeta}\right)\right).\end{aligned}$**

(iii)If $\theta\in[0,\frac{1}{2}],\operatorname*{sup}_{k\in\mathbb{N}}\frac{\gamma\eta_{u^{k-1}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k-1}}+1)+\zeta}<+\infty$ and $\operatorname*{inf}_{k\in\mathbb{N}}\frac{\zeta(\gamma-2)+1-\gamma}{2[(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k}}+1)+\zeta]}>0$, there is $k_{0}\in\mathbb{N}$ such that

• •

\begin{aligned} H_{\delta}(x^{k})-H_{\delta}(x^{\star})&=O\bigg{(}\bigg{(}\sum_{n=k_{0}}^{k-1}\frac{\gamma\eta_{u^{n}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{n}}+1)+\zeta}\bigg{)}^{\!\!\frac{-1}{1-2\theta}}\bigg{)},\end{aligned}**

• •

\begin{aligned} \|x^{\star}-x^{k}\|&=O\bigg{(}\bigg{(}\sum_{n=k_{0}}^{k-2}\frac{\gamma\eta_{u^{n}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{n}}+1)+\zeta}\bigg{)}^{\!\!\frac{-\theta}{1-2\theta}}\bigg{)}.\end{aligned}**

Proof.

The proof technique follows the route in [40]. We present the proof detail for the case $\theta\in[0,1)$, because the relation implies that the Algorithm 2 enjoys a linear convergence rate which differs from the local convergence rate analysis based on KL property $\operatorname{in}$ [19, 13, 22, 23]. Let $R_{k}:=H_{\delta_{k}}(x^{k})-H_{\delta}(x^{\star})\geq 0$, we can suppose that $R_{k}>0$ for all $k\in\mathbb{N}$, because otherwise the algorithm terminates in a finite number of steps. Since $x^{k}$ converges to $x^{\star}$, there exists $k_{0}\in\mathbb{N}$ such that, for all $k\geq k_{0}$, we have $x^{k}\in\mathcal{B}_{\rho}(\hat{x})$ where the KL inequality holds. Using successively $\mathbf{H}_{1},\mathbf{H}_{2}$ and the KL inequality, we obtain

[TABLE]

for each $k\geq k_{0}$. Let us now consider different cases for $\theta{:}$

Case $\theta=1{:}$ Suppose that $R_{k}>0$ for all $k\in\mathbb{N}.$ Then, for each $k\geq k_{0}$, we have

[TABLE]

Since $R_{k}$ converges, we must have $\inf_{k\in\mathbb{N}}\frac{[\zeta(\gamma-2)+(1-\gamma)]\gamma\eta_{u^{k}}}{2[(\Re\zeta+1)(Ls\cdot\gamma\cdot\eta_{u^{k}}+1)+\zeta]^{2}}=0$, which is a contradiction. Therefore, there exists some $k\in\mathbb{N}$ such that $R_{k}=0$, and the algorithm terminates in a finite number of steps.

Case $\theta\in[0,1]:$ $\nu_{k}=\frac{\gamma\eta_{u^{k-1}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k-1}}+1)+\zeta}$, $\bar{\nu}:=\sup_{k\in\mathbb{N}}\frac{\gamma\eta_{u^{k-1}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k-1}}+1)+\zeta}$,

$\overline{m}:=\inf_{k\in\mathbb{N}}\frac{\zeta(\gamma-2)+1-\gamma}{2[(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{k}}+1)+\zeta]}$, $c=\frac{\overline{m}}{C^{2}(1+\bar{\nu})}$ and, for each $k\in\mathbb{N},\beta_{k}:=\frac{\nu_{k}\overline{m}}{C^{2}}.$ For each $k\geq k_{0}$, (42) gives

[TABLE]

Subcase $\theta\in[\frac{1}{2},1]:$ Since $R_{k}\to 0$ and $0<2-2\theta\leq 1$, we may assume, by enlarging $k_{0}$ if necessary, that $R_{k+1}^{2-2\theta}\geq R_{k+1}$ for all $k\geq k_{0}.$ Inequality (43) implies $(R_{k}-R_{k+1})\geq\beta_{k+1}R_{k+1}$, equivalently, $R_{k+1}\leq R_{k}\left(\frac{1}{1+\beta_{k+1}}\right)$ for all $k\geq k_{0}.$ By induction, we obtain

[TABLE]

for all $k\geq k_{0}$. However, $\ln\left(\frac{1}{1+\beta_{n+1}}\right)\leq\frac{-\beta_{n+1}}{1+\beta_{n+1}}\leq\frac{-1}{1+\bar{\nu}}\beta_{n+1}$, and so

[TABLE]

Subcase $\theta\in[0,\frac{1}{2}]:$ Recall from inequality (43) that $R_{k+1}^{2\theta-2}(R_{k}-R_{k+1})\geq\beta_{k+1}.$ Setting $\phi(t):=\frac{C}{1-2\theta}{t}^{2\theta-1}$, we immediately obtain $\phi^{\prime}(t)=-Ct^{2\theta-2}$, and

[TABLE]

On the one hand, if we suppose that $R_{k+1}^{2\theta-2}\leq 2R_{k}^{2\theta-2}$, then

[TABLE]

On the other hand, we have that $R_{k+1}^{2\theta-2}>2R_{k}^{2\theta-2}$. Since $2\theta-2<2\theta-1<0$, we obtain $\frac{2\theta-1}{2\theta-2}>0$. Thus $R_{k+1}^{2\theta-1}>\Lambda R_{k}^{2\theta-1}$, where $\Lambda:=2^{\frac{2\theta-1}{2\theta-2}}>1.$ Therefore,

[TABLE]

with $C^{\prime}:=\frac{C}{1-2\theta}(\Lambda-1)R_{k_{0}}^{2\theta-1}>0.$ Since $\beta_{k+1}\leq\frac{\bar{\nu}\overline{m}}{C^{2}}$, we can write

[TABLE]

Setting $c:=\min\{\frac{C}{2},\frac{C^{\prime}C^{2}}{\bar{\nu}\overline{m}}\}>0,$ we can write $\phi(R_{k+1})-\phi(R_{k})\geq c\beta_{k+1}$ for all $k\geq k_{0}.$ This implies

[TABLE]

which is precisely $R_{k+1}\leq D\left(\sum_{n=k_{0}}^{k}\frac{\gamma\eta_{u^{n}}}{(\Re\zeta+1)(L_{s}\cdot\gamma\cdot\eta_{u^{n}}+1)+\zeta}\right)^{\frac{-1}{1-2\theta}}$ with $D=\left(\frac{c\overline{m}(1-2\theta)}{C^{3}}\right)^{\frac{-1}{1-2\theta}}.$ ∎

4 Experiments

We compare sPDOME and PDOME with widely recognized first-order methods, including PG [1] and mAPG [13], and second-order methods, such as PANOC [15] and its variant [16], as well as PDOM [17]. All algorithm comparisons were performed on a Windows 10 computer equipped with an Intel(R) Core(TM) i7-12700 2.10 GHz CPU and 16GB of memory, with all algorithms implemented and run in MATLAB. Each benchmark algorithm is fine-tuned for a fair comparison. The hyperparameter settings in the Algorithm 1 are as follows: $\eta=1/L_{s}$, $\gamma=0.98$ and $\epsilon^{\mathrm{abs}}=\epsilon^{\mathrm{rel}}=10^{-12}$. The hyperparameter settings in the Algorithm 2 are as follows: $\eta=1/L_{s}$, $\gamma=0.94$ and $\epsilon^{\mathrm{abs}}=\epsilon^{\mathrm{rel}}=10^{-12}$. All algorithms share the same randomly selected initial point $x^{0}$ and are stopped if $\left\|x^{k+1}-x^{k}\right\|/\left(1+\left\|x^{k+1}\right\|\right)<10^{-8}$ or $k>2000$. In the $k$-th iteration, we calculate the subdifferential and the normalized recovery error, NRE$(k)=\|x^{k}-x^{\star}\|/\|x^{\star}\|.$

4.1 Nonconvex Sparse Recovery Problem

The sparse recovery aims to recover the original sparse signal from an under-determined set of measurements.

[TABLE]

where $\Delta\in\mathbb{R}^{m\times n}$ denotes the measurement matrix, $\Upsilon\in\mathbb{R}^{n\times n}$ is the sparse transformation basis and $\lambda>0.$ In this context, the overall sensing matrix is specifically the subsampled Discrete Cosine Transform (DCT) [3], where $\Delta$ serves as the sub-sampling matrix and $\Upsilon$ represents the DCT basis.

The Hessian of $\bm{\mathrm{M}}$ in (44) takes the form of $\Upsilon^{\top}\Delta^{\top}\Delta\Upsilon.$ The fast matrix inversion involves expressing $\left(\Upsilon^{\top}\Delta^{\top}\Delta\Upsilon\right)^{-1}$ as $\Upsilon^{-1}\left(\Delta^{\top}\Delta\right)^{-1}{\Upsilon^{\top}}^{-1}=\Upsilon^{\top}\left(\Delta^{\top}\Delta\right)^{-1}\Upsilon$, exploiting the transpose property of the DCT basis. Notably, $\Delta^{\top}\Delta$ is a rank-deficient diagonal matrix. To ensure the invertibility, a small positive value $\iota$ is incorporated along the diagonal. Due to its diagonal structure, the computational complexity of the inversion is $\mathcal{O}(n).$

The experimental settings are summarized below: $m=$ $n/2$, ${y}={\Delta\Upsilon}x^{\star}$, $A=\Delta\Upsilon$, where the ground truth $x^{\star}$ is sparse with randomly generated entries and $\lambda=0.1\left|{A}^{T}{y}\right|_{\infty}$, following the strategy outlined in [41].

In Subfigures $(a)$ and $(b)$ of Figure 2 depict specific convergence behavior of the compared algorithm for one realization with $m=500$ in Table 2. sPDOME and PDOME outperform other baseline algorithms in terms of convergence speed and global optimality. In Subfigure $(c)$ of Figure 2, the phase transition curve is illustrated, where realizations with random initialization are considered successful if NRE$<10^{-4}.$ The results demonstrate a significantly higher success recovery rate for the sPDOME and PDOME algorithms compared to the benchmark algorithms.

4.2 Nonconvex Sparse Approximation Problem

Here, we address the challenge of identifying a sparse solution to a least-squares problem. As elaborated in [7], this is accomplished by tackling the following nonconvex optimization problem:

[TABLE]

where $\lambda>0$ is a regularization parameter, and $\|x\|_{1/2}=\left(\sum_{i=1}^{n}|x_{i}|^{1/2}\right)^{2}$ is the quasi-norm $\ell_{1/2}$, a nonconvex regularizer whose role to induce the solution of (45). Function $\|x\|_{1/2}^{1/2}$ is separable, and its proximal mapping can be computed in closed form as follows, see ([7, Theorem 1] ): for $i=1,\ldots,n$.

[TABLE]

where $p_{\eta}(x_{i})=\textbf{arccos}\left(\eta/8\left(|x_{i}|/3\right)^{-3/2}\right).$ We performed experiments using the setting of [42, Sec. 8.2]: matrix $A\in\mathbb{R}^{m\times n}$ has $m=n/5$ rows and was generated with random Gaussian entries, with zero mean and variance $1/m.$ Vector $b$ was generated as $b=Ax_{\mathrm{orig}}+v$ where $x_{\mathrm{orig}}\in\mathbb{R}^{n}$ was randomly generated with $k=5$ nonzero normally distributed entries, and $v$ is a noise vector with zero mean and variance $1/m$. Then we solved problem (45) using $x^{0}=0$ as starting iterate for all algorithms. In this experiment, Figure 3 shows that the proposed algorithms are effective.

The first two subfigures $(a)$ and $(b)$ of Figure 3 illustrate the convergence behavior of a single realization with m=500 from Table 3. Compared with the benchmark algorithms, the sPDOME and PDOME algorithms converge to the critical point faster and achieve smaller recovery errors. Subfigure $(c)$ of Figure 3 shows that the successful recovery rates of the sPDOME and PDOME algorithms are significantly higher than those of the benchmark algorithms.

Table 3 provides an overview of the average performance of the proposed algorithm and benchmark methods across various problem sizes. It is noticeable that sPDOME and PDOME converge more rapidly and exhibit a stronger capability to reach a better optimum than other algorithms. This is clearly reflected in the fact that they require fewer iterations to get close to the critical point and achieve a smaller NER.

5 Conclusion

In this paper, the PDOME algorithm are proposed for nonconvex and nonsmooth problems with a quadratic term. During the iteration process of PDOME algorithm, constructs and minimizes a majorant function in a hybrid direction based on extrapolation. Theoretical analysis confirms that the algorithm can achieve convergence to a critical point, and its global convergence rate is studied based on the KL property. Numerical experiments show that the sPDOME and PDOME algorithms converges faster in nonconvex problems and can better converge to a local optimal solution. Building on the research in this paper, future work will further expand the content related to quadratic functions and Bregman distances at the algorithm level.

Declarations

Funding: This research is supported by the National Natural Science Foundation of China (NSFC) grants 92473208, 12401415, the Key Program of National Natural Science of China 12331011, the 111 Project (No. D23017), the Natural Science Foundation of Hunan Province (No. 2025JJ60009), the Tianchi Talent Program of Xinjiang Uygur Autonomous Region (CZ001328).

Data Availability: Enquiries about data/code availability should be directed to the authors.

Competing interests: The authors have no competing interests to declare that are relevant to the content of this paper.