Paper deep dive
Maximum Tsallis Entropy Distributions for Robust and Efficient Sparse Learning from Correlated Data
Kai Yang, Masoud Asgharian, Celia M. T. Greenwood
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/19/2026, 4:29:58 AM
Summary
This paper proposes using the qGaussian distribution, derived from maximizing Tsallis entropy, as a robust alternative to Gaussian distributions for statistical sparse learning, particularly for correlated and heterogeneous data in biostatistics. The authors re-derive the multivariate qGaussian probability density function and introduce a novel optimization framework that adapts numerical methods for finding equilibria in flows to solve composite optimization problems. Specifically, they apply this framework to the Hager-Zhang conjugate gradient algorithm to develop a numerically stable and efficient proximal conjugate gradient algorithm for sparse statistical learning.
Entities (8)
Relation Signals (6)
qGaussian distribution → derivedfrom → Tsallis entropy
confidence 98% · we propose the use of the qGaussian distribution, derived from Tsallis entropy maximization
Hager-Zhang conjugate gradient algorithm → adaptedfor → sparse statistical learning
confidence 96% · Applying this framework to the Hager-Zhang conjugate gradient algorithm... we develop a numerically stable and efficient algorithm for sparse statistical learning.
qGaussian distribution → alternativeto → Gaussian distribution
confidence 95% · presents a viable and flexible alternative to Gaussian-based methods
qGaussian distribution → appliedin → biostatistics
confidence 92% · This is notably relevant in biostatistics... we propose the use of the qGaussian distribution
Tsallis entropy → generalizes → Shannon's entropy
confidence 90% · constructed an entropy similar to Shannon’s entropy... Tsallis entropy is also known as non-extensive entropy
Moreau envelope → usedin → optimization framework
confidence 88% · By employing the Moreau envelope and linearizing the smooth term, we develop a framework
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:This paper addresses the limitations of Gaussian distribution assumptions in statistical sparse learning, particularly in modeling correlated and heterogeneous data. Conventional Gaussian models often lack robustness towards outliers and underlying distribution assumptions. To overcome these limitations, we propose the use of the $q$Gaussian distribution, derived from Tsallis entropy maximization, as a robust alternative. This is notably relevant in biostatistics, where the presence of correlated observations and heterogeneity, such as in genetic and longitudinal studies, is prevalent. Our contributions include modeling of correlated data through the re-derived multivariate probability density function from Tsallis entropy maximization, thereby addressing the limitations inherent in conventional Gaussian models. Furthermore, we introduce a novel framework that adapts numerical methods designed to find equilibria in flows to tackle composite optimization problems prevalent in statistical sparse learning. Applying this framework to the Hager-Zhang conjugate gradient algorithm \cite{Hager2005}, we develop a numerically stable and efficient algorithm for sparse statistical learning. The $q$Gaussian distribution, informed by the principle of maximizing Tsallis entropy, presents a viable and flexible alternative to Gaussian-based methods. This paper not only contributes to the theoretical understanding of statistical distributions and optimization techniques, but also paves the way for practical data analysis.
Tags
Links
- Source: https://arxiv.org/abs/2608.17244v1
- Canonical: https://arxiv.org/abs/2608.17244v1
Trouble viewing inline? Open PDF directly →
Full Text
144,709 characters extracted from source content.
Expand or collapse full text
Maximum Tsallis Entropy Distributions for Robust and Efficient Sparse Learning from Correlated Data ang kai.yang2@mail.mcgill.ca ud Asgharian masoud.asgharian2@mcgill.ca a M.T. Greenwood celia.greenwood@mcgill.ca Abstract This paper addresses the limitations of Gaussian distribution assumptions in statistical sparse learning, particularly in modeling correlated and heterogeneous data. Conventional Gaussian models often lack robustness towards outliers and underlying distribution assumptions. To overcome these limitations, we propose the use of the qqGaussian distribution, derived from Tsallis entropy maximization, as a robust alternative. This is notably relevant in biostatistics, where the presence of correlated observations and heterogeneity, such as in genetic and longitudinal studies, is prevalent. Our contributions include modeling of correlated data through the re-derived multivariate probability density function from Tsallis entropy maximization, thereby addressing the limitations inherent in conventional Gaussian models. Furthermore, we introduce a novel framework that adapts numerical methods designed to find equilibria in flows to tackle composite optimization problems prevalent in statistical sparse learning. Applying this framework to the Hager-Zhang conjugate gradient algorithm 22, we develop a numerically stable and efficient algorithm for sparse statistical learning. The qqGaussian distribution, informed by the principle of maximizing Tsallis entropy, presents a viable and flexible alternative to Gaussian-based methods. This paper not only contributes to the theoretical understanding of statistical distributions and optimization techniques, but also paves the way for practical data analysis. †shortheadings: MTE Distributions for Robust and Efficient Sparse Learning / Yang, Asgharian, and Greenwood†firstpage: 1†editor: My editor keywords Statistical Learning, high–dimensional, Entropy, Robust, Optimization 1 Introduction In the realm of statistical sparse learning, the pursuit of robust and efficient methodologies remains paramount, especially when confronted with the complexities of correlated data. The principle of maximizing Shannon’s entropy stands as a pivotal framework that has led to the derivation of nearly all frequently utilized statistical distributions to date 12. This principle’s application has notably revealed that the multivariate Gaussian distribution maximizes Shannon’s entropy under the first two moments, significantly influencing the landscape of statistical sparse learning. The Gaussian assumption has become a fundamental cornerstone of numerous statistical sparse learning problem formulations, and its core assumptions are rarely re-examined or challenged. However, the Gaussian distribution’s features, particularly its exponential tail decay and the absence of a shape parameter, can present substantial limitations. Specifically, the lack of robustness towards outliers and a limited capacity to accurately represent the distribution’s shape results in violations of the Gaussian assumption in statistical modeling. Such violations have many practical repercussions, including the potential for erroneous Type I error rates and the lack of robustness towards distribution shape when estimating the dispersion parameter, motivating the development of alternative approaches. Dispersion or volatility parameters encapsulate critical and often decisive information about distributions. Their estimations, specifically in transformations of predicted outcomes, are often indispensable. For instance, in the context of log-normal distributions, the mean is directly influenced by the volatility parameter derived from the underlying Gaussian distribution. Likewise, principles like the Law of the Unconscious Statistician (LOTUS), which rely on accurate estimation of volatility and precise understanding of the distribution’s shape, highlight the importance of determining this parameter for dependable prediction and statistical modeling. Volatility estimation is of great importance in finance. Specifically, Itô’s lemma, often used in stochastic calculus for option pricing, explicitly highlights the importance of volatility’s contribution. Delving into the realm of stochastic calculus, Itô’s lemma provides a mathematical framework that elegantly captures volatility’s impact on dynamic systems. Itô’s lemma states that for a twice-differentiable function f, df(t,Xt)=(∂f∂t+μ∂f∂x+12σ2∂2f∂x2)dt+σ∂f∂xdWt,df (t,X_t )= ( ∂ f∂ t+μ ∂ f∂ x+ 12σ^2 ∂^2f∂ x^2 )dt+σ ∂ f∂ xdW_t, (1) where the term σ2∂2f∂x2σ^2 ∂^2f∂ x^2 specifically denotes the contribution of volatility to changes in the function f. This mathematical representation is pivotal in finance, where the phenomenon, termed volatility smile, challenges the foundational assumptions of the Black–Scholes–Merton model, signaling empirical deviations from expected normality in option pricing models. These deviations have propelled the exploration of alternative distributions capable of more accurately reflecting market realities 42. In response to these limitations of Gaussian distributions, the qqGaussian distribution, derived from the maximization of the Tsallis entropy, emerges as a compelling alternative. The qqGaussian distribution is celebrated for its flexibility in modeling the diverse shapes of bell-curved distributions, including the ability to account for heavy-tailed distributions — a feature crucial for robust modeling of financial returns. It provides a more accurate representation of financial returns on platforms such as the New York Stock Exchange and National Association of Securities Dealers Automated Quotations 5; 6; 15. Despite its proven advantages in finance, the incorporation of Tsallis entropy-maximizing distributions within the domain of statistical sparse learning and biostatistics remains limited. To the best of our knowledge, this paper represents the initial endeavor to apply Tsallis entropy-maximizing distributions for biostatistical data modeling. Correlated observations, frequently encountered in genetic and longitudinal studies 18; 48; 14, as well as heterogeneity of the variance, will be specifically addressed by our proposed Tsallis entropy-maximizing model for correlated data. Hence, this paper advocates for the application of the qqGaussian distribution in modeling correlated data within sparse statistical learning frameworks. Our approach relaxes the conventional reliance on normality assumptions; We aim to demonstrate that the intricate characteristics of qqGaussian distributions can profoundly enhance the modeling of correlated data, offering a robust and versatile alternative to conventional Gaussian-based methods. Maximum likelihood estimation is one of the most commonly used estimation techniques. However, the estimation process encounters notable computational obstacles when dealing with high–dimensional and extra-large datasets. Oracle penalties, favored for their efficacy in facilitating variable selection, present an attractive yet complex solution to sparse learning problems. However, oracle penalties are notable for their nonconvex and nonsmooth nature 40, which lead to considerable optimization challenges. Recently, proximal methods have demonstrated an unmatched speed of convergence, thereby surpassing most other approaches in efficiently handling estimation in nonsmooth problems 26. Simultaneously, the Krylov subspace method, recognized among the top ten algorithms for computing in science and engineering of the twentieth century, lays a solid foundation for numerical analysis. The conjugate gradient method, a prominent member of the Krylov subspace methods, has been applied extensively in various areas and is a fundamental numerical tool in solving partial differential equations 41. While the nonlinear conjugate gradient performs exceptionally well in terms of its convergence speed and numerical stability, much better than accelerated gradient or gradient descent, its global convergence depends on the line–search step, whereas accelerated gradient and gradient descent methods do not necessarily require the line search step to achieve global convergence 20; 63. Motivated by these methodologies, our paper introduces a proximal conjugate gradient method that can be applied to solve qqGaussian sparse learning problems. This method aims to effectively combine the theoretical strengths of both proximal methods and Krylov subspace techniques. Additionally, our paper addresses the line search step needed for the proximal nonlinear conjugate gradient method to establish global convergence. We re-derive the probability density function for the multivariate qqGaussian distribution from a Tsallis entropy maximizing perspective in Lemma 1. This derivation allows for a nuanced understanding and application of this model in statistical analysis. Furthermore, our contributions in this paper are as following: 1. We apply the derived density to model correlated and heterogeneous data effectively, while carrying out the sparse statistical learning at the same time. 2. Sparse statistical learning involves minimizing a composite optimization problem, aimed at minimizing a composite objective function composed of a globally Lipschitz-smooth term, which may be nonconvex, and a convex nonsmooth term. A variety of numerical methods are available to find equilibrium points for globally Lipschitz flows. By employing the Moreau envelope and linearizing the smooth term, we develop a framework that allows any numerical method designed for finding equilibrium points in globally Lipschitz flows to be adapted into a numerical optimization algorithm for minimizing the composite objective function. 3. Leveraging the framework introduced above, we implement it with the state-of-the-art Hager-Zhang conjugate gradient method 22. This implementation yields a proximal conjugate gradient algorithm that is not only computationally efficient but also numerically stable, suitable for a wide range of statistical sparse learning challenges. This includes the robust sparse learning approach we devised based on the concept of maximizing the Tsallis entropy distribution. The structure of the paper is organized as follows: Section 2 delves into the foundational properties of Tsallis entropy, drawing upon previous literature to establish a comprehensive background. Following this, Section 3 introduces the concept of q−q-moments. This section then elaborates on employing Tsallis entropy maximizing distribution to effectively model the q−q-correlation structure. In Section 3, we also re-derive the probability density function maximizing Tsallis entropy under the first and second central q−q-moment constraints, incorporating all relevant parameters for a likelihood-based approach to statistical analysis. Our discussion transitions to the challenges and strategies of optimization in Section 4. This section is twofold; initially, in Section 4.1, we present essential background knowledge from variational and nonsmooth analysis. This foundation is critical for our novel contribution: the development of a proximal framework to transform any first-order numerical optimization algorithm to a proximal counterpart by leveraging the properties of the Moreau envelope, detailed in Section 4.2. In Section 4.3, we apply this innovative framework to the state-of-the-art Hager-Zhang conjugate gradient algorithm. This adaptation produces a proximal version for tackling sparse statistical learning challenges. The efficacy of this method is further showcased in Section 5, where we outline the application of our proximal Hager-Zhang conjugate gradient algorithm to optimize a penalized qqGaussian likelihood function. This section also lays out a map from problem formulation to the practical aspects of prediction using models trained with our approach. Finally, Section 6 synthesizes our contributions, offering a reflective conclusion and proposing avenues for future research. 2 Tsallis Entropy For an arbitrary random variable X, Shannon’s Entropy 50 poses the definition H(X)≔−log(p(X))=−∫log(p(x))dμX,H (X ) -E (p (X ) )=- (p (x ) )d _X, (2) where p is the likelihood function for X. Over a given (likelihood) function space ≔p(x)|∀x∈,p(x)≥0,|log(p(X))|<∞ and ∫1dμX=1,P \p (x )|∀ x ,\ p (x )≥ 0, |E (p (X ) ) |<∞ and _X1d _X=1 \, the Shannon’s entropy is a strictly concave function, which implies uniqueness of the maximizer. Many commonly-used distributions have been shown to maximize Shannon’s entropy under certain given constraints 12. For example, uniform distribution, whether in discrete or continuous case , are to maximize (2) over a compact support, with open sets defined by discrete or Euclidean topology, respectively. The exponential distribution is defined as maximizing (2) over ℝ≥0R_≥ 0 and with a constraint that the first moment is a constant, 1λ 1λ; where λ later turns out to be the scale parameter. And the Gaussian distribution maximizes (2) over ℝR with given mean and variance. More examples can be given. For example, the constraint to obtain a Laplace distribution is a given mean absolute deviation, etc. Additivity is a key element of Shannon’s entropy. That is, let A1,A2A_1,A_2 be two independent event sets, then the information of the intersection I(ℙ(A1∩A2))=I(ℙ(A1)⋅ℙ(A2))=I(ℙ(A1))+I(ℙ(A2))I (P (A_1∩ A_2 ) )=I (P (A_1 )·P (A_2 ) )=I (P (A_1 ) )+I (P (A_2 ) ) — such homomorphism was considered particularly useful in Shannon’s view 50. Later in the 1980s, 56 constructed an entropy similar to Shannon’s entropy but without the additivity property. To see how Tsallis’ entropy was developed, first we look at Tsallis’ q−q-exponential function, which is defined as expq:ℝ↦ℝ _q:R , given by expqx≔((1+(1−q)x))11−q;1+(1−q)x>00;else _qx cases ( (1+ (1-q )x ) ) 11-q;&1+ (1-q )x>0\\ 0;&else cases (3) For q>1q>1, expq _q is bijective over (0,1q−1) (0, 1q-1 ). The inverse function, called the q−q-logarithmic function, is given by lnqx≔x1−q−11−q. _qx x^1-q-11-q. (4) Based on this deformed q−q-exponential function, 56 developed Tsallis entropy by replacing the log function in 2 with q−logq- function (4) and replacing the expectation with q−q-expectation operator 56: Sq(X) S_q (X ) =−∫pq(x)lnqp(x)dx≕−qlnqp(X) =- _Xp^q (x ) _qp (x )dx -E_q _qp (X ) (5) =1q−1(1−∫pq(x)x) = 1q-1 (1- _Xp^q (x )dx ) (6) where q∈ℝ∖1q \1\ is a constant, and qf(X)≔∫f(x)⋅pq(x)x=⟨f(x),(dμX)q⟩E_qf (X ) _Xf (x )· p^q (x )dx= f (x ), (d _X )^q is referred to as the q−q-expectation operator. Tsallis entropy is also known as non-extensive entropy; namely for arbitrary independent two random variable X1,X2X_1,X_2: Sq(X1,X2)=Sq(X1)⊕qSq(X2),S_q(X_1,X_2)=S_q(X_1) _qS_q(X_2), (7) where “⊕q _q” is defined as ∀a,b∈ℝ∀ a,b , a⊕qb≔a+b+(1−q)ab.a _qb a+b+ (1-q )ab. (8) Expectation has been used to characterize statistical distributions. However, one significant drawback of the expectation (linear) operator is the lack of continuity for some distributions; such as the Cauchy distribution. Therefore, the q−q-expectation operator, qE_q, provides robustness when characterizing the distributions in the real domain. If the tail of the function vanishes at a rate of O((logx)−1)O ( ( x )^-1 ), the function will not have a proper integral if the support is unbounded. Thus, for any distribution whose likelihood function is bounded in uniform norm, ∃q∈ℝ>0∃ q _>0 such that qE_q is continuous at the likelihood function in the function space we are considering. 3 Tsallis Entropy Maximizing Distribution to Accommodate the q−q-Correlation Structure The Gaussian distribution maximizes the Shannon’s entropy in the following problem: maxϕ∈ \ _φ −∫ℝnϕ(x)log(ϕ(x))dx - _R^nφ (x ) (φ (x ) )dx s.t. ϕ≥0; φ≥ 0; ∫ℝnϕ(x)x=1; _R^nφ (x )dx=1; ∫ℝnx⋅ϕ(x)x=0; _R^nx·φ (x )dx=0; (9) ∫ℝnxxT⋅ϕ(x)x=Σ; _R^nx^T·φ (x )dx= ; for some n∈ℕ+n _+ and Σ∈ℝn×n,Σ≻0 ^n× n,\ 0. For the sake of parsimony, in (9) we assume that the distribution is centered. To set the central trend parameter, or the mean parameter in the specific case of the Gaussian distribution, the likelihood function ϕφ can be simply translated x↦x−μx x-μ to incorporate the parameter μ for the central trend. Entropy functions are invariant under translation. Similarly to how the multivariate Gaussian distribution maximizes Shannon’s entropy in a Euclidean space, the multivariate qqGaussian distribution maximizes Tsallis entropy in a Euclidean space. Specifically, the optimization problem is formulated as: maxϕ∈Lq(ℝn) \ _φ∈ L^q (R^n ) −∫ℝnϕq(x)dx - _R^nφ^q (x )dx (10) s.t. ϕ≥0; φ≥ 0; ∫ℝnϕ(x)x=1; _R^nφ (x )dx=1; (11) ∫ℝnx⋅ϕq(x)x∫ℝnϕq(x)x=0; _R^nx·φ^q (x )dx _R^nφ^q (x )dx=0; (12) ∫ℝnxxT⋅ϕq(x)x∫ℝnϕq(x)x=Σ, _R^nx^T·φ^q (x )dx _R^nφ^q (x )dx= , (13) where q>1q>1. The feasible set of Lebesgue space Lq(ℝn)L^q (R^n ) is to ensure the well–definedness of Tsallis entropy. The normalization constraint (11) implies that ϕ∈L1φ∈ L^1; however, ϕ∈L1φ∈ L^1 does not imply ϕ∈Lqφ∈ L^q, as the embedding property L1⊆LqL^1 L^q fails to hold for Lebesgue measure on ℝnR^n. As an example, consider the one-dimensional example of the probability density function ϕ~(x)=14|x|−12for x∈(−1,1)∖0;0else. φ (x )= cases 14 |x |^- 12&for x∈ (-1,1 ) \0 \;\\ 0&else. cases (14) Clearly, ϕ~∈L1 φ∈ L^1 but ϕ~∉L2 φ ∈ L^2. Note that (12) and (13) are the first and second moment constraints using the q−q-expectation operator qE_q. As noted by 34, maximizing any member of the generalized class of power-law entropies, including Renyi entropy, Havrda and Charvat entropy, Arimoto entropy, and Tsallis entropy, all yield the identical power-law objective function (10). Regarding the constraints, ∫ℝnϕq(x)x _R^nφ^q (x )dx is the normalization factor for the q−q-expectation. 34 further noted that optimizing the problem formulated above is equivalent to the following problem: maxφ∈Ls(ℝn) \ _ ∈ L^s (R^n ) ∫ℝnφs(x)x _R^n ^s (x )dx (15) s.t. φ≥0; ≥ 0; ∫ℝnφ(x)x=1; _R^n (x )dx=1; ∫ℝnx⋅φ(x)x=0; _R^nx· (x )dx=0; ∫ℝnxxT⋅φ(x)x=Σ. _R^nx^T· (x )dx= . In (15), s≔q−1∈(0,1)s q^-1∈ (0,1 ), thus Ls(ℝn)L^s (R^n ) is a quasi-normed space; φ(x)≔ϕq(x)∫ℝnϕq(x)x. (x ) φ^q (x ) _R^nφ^q (x )dx. If the maximizer of (15) is φ , then the maximizer of (10), ϕφ, will be normalized ϕ(x)∝φ1/q.φ (x ) ^1/q. (16) Several important properties were proposed previously regarding the qqGaussian distributions in previous studies 57; 11. Notably, 1. Using Bregman information divergence, Problem (15) has a unique maximizer of the form φ(x,s)=As(1−(s−1)β′⟨x,Σ−1x⟩)+1s−1 (x;s )=A_s (1- (s-1 )β x, ^-1x )_+ 1s-1 (17) for some s∈(n+2,∞)∖1s∈ ( nn+2,∞ ) \1 \, normalization constant ArA_r, and some dispersion parameter β′β . 2. If X∼qGaussian(q,Σ)X qGaussian (q, ), H∈ℝn~×nH n× n and rank(H)=n~rank (H )= n. Then X~∼qGaussian(q~,HΣHT) X qGaussian ( q,H H^T ) with 21−q~−1−n~=21−q−1−n. 21- q^-1- n= 21-q^-1-n. (18) 3. If X1,X2X_1,X_2 are both qqGaussian random vectors but independent, a linear combination of H1X1+H2X2H_1X_1+H_2X_2 is not qqGaussian. 4. The duality property: if X∼qGaussian(q,Σ)X qGaussian (q, ) with 1<q<1+2n1<q<1+ 2n, let the degree of freedom for X be m≔2q−1−nm 2q-1-n and Λ≔mΣ m , then X1−⟨X,Λ−1X⟩ X 1- X, ^-1X ∼qGaussian(q~,m+4Σ) qGaussian ( q, mm+4 ) with 1q~−1−1 1 q^-1-1 =11−q−1−n2−1, = 11-q^-1- n2-1, and 0<q~<10< q<1. Property 1 will be used in our Lemma 1. Property 2 implies that any components of a qqGaussian random vector are also qqGaussian, while Property 3 implies that two independent qqGaussian vectors are not jointly qqGaussian. By the equivalence of problems (10) and (15) discussed before, (17) can be rewritten as ϕ(x,q,Σ)=(α−β⟨x,Σ−1x⟩)+11−qφ (x;q, )= (α-β x, ^-1x )_+ 11-q (19) for some constant (parameter) α,βα,β, q∈(0,1+2n)∖1q∈ (0,1+ 2n ) \1 \; x+≔max(0,x)x_+ \ (0,x ). As shown later in the proof of Lemma 1, the dimension-related upper bound 1+2n1+ 2n is due to the normalization constraint (11). When 0<q<10<q<1, the density represents a distribution with bounded support; when q>1q>1, the density is a generalization of the bell curve distributions, and with q↘1q 1 the Gaussian distribution is recovered. A higher value of q corresponds to heavier tails in shape. The duality between the qqGaussian random vectors with 0<q<10<q<1 and 1<q<1+2n1<q<1+ 2n was given by 57, which we discussed in Property 4 in Section 3. Distributions with bounded support correspond to 0<q<10<q<1, and distributions with heavy tails correspond to 1<q<1+2n1<q<1+ 2n. For the scope of this paper, we will focus only on the heavy-tail distributions; i.e., the case when q>1q>1. 60 derived the qqGaussian probability density function for 1<q<n+4n+21<q< n+4n+2, when the multivariate qqGaussian density becomes the scaled density of the multivariate student’s t distribution. To incorporate the case of q∈[1+2n+2,1+2n)q∈[1+ 2n+2,1+ 2n), when the variance does not exist but the q−q-variance can be used to capture the volatility/dispersion of the data, 59 also derived the resulting density; however, since a typo was found in that paper, we re-derive the density in Lemma 1. The parameters presented in the density formula (20) are of particular interest to statisticians, as parameter inference is the key to statistical analysis and prediction. The case of q∈[1+2n+2,1+2n)q∈[1+ 2n+2,1+ 2n) will allow the resulting qqGaussian distribution to incorporate the wider class of distributions without finite moments but finite q−q-moments; such as the Cauchy distribution. Therefore, modeling using the qqGaussian distribution with q allowed to take the value in [1+2n+2,1+2n)[1+ 2n+2,1+ 2n) will be more robust. Lemma 1. When q∈(1,1+2n)q∈ (1,1+ 2n ), the unique solution to (10) is: p(x,q,Σ)=1|πΣ|1/2⋅Γ(1q−1)Γ(1q−1−n2)⋅(2q−1−n)−n2⋅(1+(2q−1−n)−1⋅⟨x,Σ−1x⟩)11−q.p (x;q, )= 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 )· ( 2q-1-n )^- n2· (1+ ( 2q-1-n )^-1· x, ^-1x ) 11-q. (20) proof By (17) and the equivalence of the problems (10) and (15), let the solution to (10) be denoted by p(x,q,Σ)=1Z(γ+⟨x,Σ−1x⟩)11−qp (x;q, )= 1Z (γ+ x, ^-1x ) 11-q (21) for some Z,γ>0Z,γ>0. Feasibility for problem 10 when q∈(1,1+2n)q∈ (1,1+ 2n ) was given in 57. Hence, the strictly concavity of the objective function 10 implies that the optimal solution is unique. The symmetry of p(x,q,Σ)p (x;q, ) is implied by (21); thus, we reformulate the problem (10) as the following equivalent problem: maxp∈Lq(ℝn) \ _p∈ L^q (R^n ) −∫ℝn(p(x;q,Σ))qdx - _R^n (p (x;q, ) )^qdx s.t. p(x,q,Σ)≥0; p (x;q, )≥ 0; ∫ℝ>0np(x,q,Σ)x=2−n; _R_>0^np (x;q, )dx=2^-n; ∫ℝnx⋅pq(x,q,Σ)x∫ℝnpq(x,q,Σ)x=0; _R^nx· p^q (x;q, )dx _R^np^q (x;q, )dx=0; ∫ℝnxxT⋅pq(x,q,Σ)x∫ℝnpq(x,q,Σ)x=Σ. _R^nx^T· p^q (x;q, )dx _R^np^q (x;q, )dx= . (22) Thus, Z Z =2n∫ℝ>0n(γ+⟨x,Σ−1x⟩)11−qx =2^n _R_>0^n (γ+ x, ^-1x ) 11-qdx =2n|Σ1/2|∫ℝ>0n(γ+⟨x,x⟩)11−qx =2^n | 12 | _R_>0^n (γ+ x,x ) 11-qdx =2n|Σ|1/2∫0∞rn−1(γ+r2)11−q(∏i=1n−2∫0π2sinn−1−i(θi)dθi⋅∫0π21θ)r =2^n | | 12 _0^∞r^n-1 (γ+r^2 ) 11-q ( _i=1^n-2 _0 π2 ^n-1-i( _i)d _i· _0 π21dθ )dr =2n|Σ|1/2⋅(∏i=1n−2∫0π2sinn−1−i(θi)dθi⋅∫0π21θ)⋅∫0∞rn−1(γ+r2)11−qr =2^n | | 12· ( _i=1^n-2 _0 π2 ^n-1-i( _i)d _i· _0 π21dθ )· _0^∞r^n-1 (γ+r^2 ) 11-qdr =π2n−1|Σ|1/2⋅(∏i=1n−212Γ(n−i2)πΓ(n−i+12))⋅∫0∞rn−1(γ+r2)11−qr =π 2^n-1 | | 12· ( _i=1^n-2 12 ( n-i2 ) π ( n-i+12 ) )· _0^∞r^n-1 (γ+r^2 ) 11-qdr =2πn2|Σ|1/2⋅(Γ(n2))−1⋅∫0∞rn−1(γ+r2)11−qr =2π n2 | | 12· ( ( n2 ) )^-1· _0^∞r^n-1 (γ+r^2 ) 11-qdr =2πn2|Σ|1/2⋅(Γ(n2))−1⋅∫0∞γn−12+11−q(rγ)n−1(1+(rγ)2)11−qr =2π n2 | | 12· ( ( n2 ) )^-1· _0^∞γ n-12+ 11-q ( r γ )^n-1 (1+ ( r γ )^2 ) 11-qdr =2πn2|Σ|1/2⋅(Γ(n2))−1⋅γn2+11−q⋅∫0∞(r′)n−1(1+(r′)2)11−qdr′ =2π n2 | | 12· ( ( n2 ) )^-1·γ n2+ 11-q· _0^∞ (r )^n-1 (1+ (r )^2 ) 11-qdr =2πn2|Σ|1/2⋅(Γ(n2))−1⋅γn2+11−q⋅∫0∞((r′)1−n(1+(r′)2)1q−1)−1dr′ =2π n2 | | 12· ( ( n2 ) )^-1·γ n2+ 11-q· _0^∞ ( (r )^1-n (1+ (r )^2 ) 1q-1 )^-1dr =2πn2|Σ|1/2⋅(Γ(n2))−1⋅γn2+11−q⋅12B(1q−1−n2,n2) =2π n2 | | 12· ( ( n2 ) )^-1·γ n2+ 11-q· 12B ( 1q-1- n2, n2 ) (23) =πn2|Σ|1/2⋅(Γ(n2))−1⋅γn2+11−q⋅Γ(1q−1−n2)Γ(n2)Γ(1q−1) =π n2 | | 12· ( ( n2 ) )^-1·γ n2+ 11-q· ( 1q-1- n2 ) ( n2 ) ( 1q-1 ) =πn2|Σ|1/2Γ(1q−1−n2)Γ(1q−1)⋅γn2+11−q. =π n2 | | 12 ( 1q-1- n2 ) ( 1q-1 )·γ n2+ 11-q. In step (23), we use the following formula for Beta function 34: ∫0∞(xα(1+xλ)β)−1x=1λB(β−1−αλ,1−αλ), _0^∞ (x^α (1+x^λ )^β )^-1dx= 1λB (β- 1-αλ, 1-αλ ), (24) where α<1,λ>0,β>0,λβ>1−α<1,\ λ>0,\ β>0,\ λβ>1-α. Well–definedness of Z and (23) implies that 1q−1−n2>0 1q-1- n2>0; i.e., q<1+2nq<1+ 2n, which is the reason for the upper bound for the choice of q. (22) implies that tr(Σ−1∫ℝnxxT⋅pq(x)x∫ℝnpq(x)x)=tr(Σ−1Σ)=n.tr ( ^-1 _R^nx^T· p^q (x )dx _R^np^q (x )dx )=tr ( ^-1 )=n. (25) Hence, since p(x)p (x ) is symmetric, tr(Σ−1∫ℝnxxT⋅pq(x,q,Σ)x∫ℝnpq(x,q,Σ)x) ( ^-1 _R^nx^T· p^q (x;q, )dx _R^np^q (x;q, )dx ) (26) =tr(Σ−1∫ℝ>0nxxT⋅pq(x,q,Σ)x∫ℝ>0npq(x,q,Σ)x) =tr ( ^-1 _R_>0^nx^T· p^q (x;q, )dx _R_>0^np^q (x;q, )dx ) =tr(∫ℝ>0nΣ−1xxT⋅pq(x,q,Σ)x)∫ℝ>0npq(x,q,Σ)x = tr ( _R_>0^n ^-1x^T· p^q (x;q, )dx ) _R_>0^np^q (x;q, )dx =∫ℝ>0ntr(Σ−1xxT⋅pq(x,q,Σ))x∫ℝ>0npq(x,q,Σ)x = _R_>0^ntr ( ^-1x^T· p^q (x;q, ) )dx _R_>0^np^q (x;q, )dx =∫ℝ>0ntr(xTΣ−1x)⋅pq(x,q,Σ)x∫ℝ>0npq(x,q,Σ)x = _R_>0^ntr (x^T ^-1x )· p^q (x;q, )dx _R_>0^np^q (x;q, )dx =∫ℝ>0n⟨x,Σ−1x⟩⋅(1Z(γ+⟨x,Σ−1x⟩)11−q)qx∫ℝ>0n(1Z(γ+⟨x,Σ−1x⟩)11−q)qx = _R_>0^n x, ^-1x · ( 1Z (γ+ x, ^-1x ) 11-q )^qdx _R_>0^n ( 1Z (γ+ x, ^-1x ) 11-q )^qdx =∫ℝ>0n⟨x,Σ−1x⟩⋅(γ+⟨x,Σ−1x⟩)q1−qx∫ℝ>0n(γ+⟨x,Σ−1x⟩)q1−qx = _R_>0^n x, ^-1x · (γ+ x, ^-1x ) q1-qdx _R_>0^n (γ+ x, ^-1x ) q1-qdx =∫ℝ>0n|Σ|1/2⟨x,x⟩⋅(γ+⟨x,x⟩)q1−qx∫ℝ>0n|Σ|1/2(γ+⟨x,x⟩)q1−qx = _R_>0^n | | 12 x,x · (γ+ x,x ) q1-qdx _R_>0^n | | 12 (γ+ x,x ) q1-qdx =∫ℝ>0n⟨x,x⟩⋅(γ+⟨x,x⟩)q1−qx∫ℝ>0n(γ+⟨x,x⟩)q1−qx = _R_>0^n x,x · (γ+ x,x ) q1-qdx _R_>0^n (γ+ x,x ) q1-qdx =∫0∞rn−1⋅r2⋅(γ+r2)q1−q⋅(∏i=1n−2∫0π2sinn−1−i(θi)dθi⋅∫0π21θ)r∫0∞rn−1⋅(γ+r2)q1−q⋅(∏i=1n−2∫0π2sinn−1−i(θi)dθi⋅∫0π21θ)r = _0^∞r^n-1· r^2· (γ+r^2 ) q1-q· ( _i=1^n-2 _0 π2 ^n-1-i( _i)d _i· _0 π21dθ )dr _0^∞r^n-1· (γ+r^2 ) q1-q· ( _i=1^n-2 _0 π2 ^n-1-i( _i)d _i· _0 π21dθ )dr =∫0∞rn+1⋅(γ+r2)q1−qr∫0∞rn−1⋅(γ+r2)q1−qr = _0^∞r^n+1· (γ+r^2 ) q1-qdr _0^∞r^n-1· (γ+r^2 ) q1-qdr =γ∫0∞((r′)−n−1⋅(1+(r′)2)q−1)−1dr′∫0∞((r′)1−n⋅(1+(r′)2)q−1)−1dr′ = γ _0^∞ ( (r )^-n-1· (1+ (r )^2 ) qq-1 )^-1dr _0^∞ ( (r )^1-n· (1+ (r )^2 ) qq-1 )^-1dr =γB(q−1−n+22,n+22)B(q−1−n2,n2) = γ B ( qq-1- n+22, n+22 )B ( qq-1- n2, n2 ) (27) =γ⋅Γ(q−1−n2−1)Γ(n+22)Γ(q−1)/Γ(q−1−n2)Γ(n2)Γ(q−1) =γ· ( qq-1- n2-1 ) ( n+22 ) ( qq-1 )/ ( qq-1- n2 ) ( n2 ) ( qq-1 ) =γ⋅n2⋅(q−1−n2−1)−1. =γ· n2· ( qq-1- n2-1 )^-1. (28) In step (27), we used (24). Combining (25) and (28), we have γ⋅n2⋅(q−1−n2−1)−1=n,γ· n2· ( qq-1- n2-1 )^-1=n, (29) which gives that γ=2q−1−n−2=2q−1−n.γ= 2qq-1-n-2= 2q-1-n. (30) Thus, the probability density function that maximizes problem (10) is: p(x,q,Σ) p (x;q, ) =(πn2|Σ|1/2Γ(1q−1−n2)Γ(1q−1)⋅(2q−1−n)n2+11−q)−1((2q−1−n)+⟨x,Σ−1x⟩)11−q = (π n2 | | 12 ( 1q-1- n2 ) ( 1q-1 )· ( 2q-1-n ) n2+ 11-q )^-1 ( ( 2q-1-n )+ x, ^-1x ) 11-q =(πn2|Σ|1/2Γ(1q−1−n2)Γ(1q−1)⋅(2q−1−n)n2)−1(1+(2q−1−n)−1⟨x,Σ−1x⟩)11−q = (π n2 | | 12 ( 1q-1- n2 ) ( 1q-1 )· ( 2q-1-n ) n2 )^-1 (1+ ( 2q-1-n )^-1 x, ^-1x ) 11-q =1|πΣ|1/2⋅Γ(1q−1)Γ(1q−1−n2)⋅(2q−1−n)−n2⋅(1+(2q−1−n)−1⋅⟨x,Σ−1x⟩)11−q. = 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 )· ( 2q-1-n )^- n2· (1+ ( 2q-1-n )^-1· x, ^-1x ) 11-q. When q≥1+2nq≥ 1+ 2n, the solution to problem (10) does not exist, due to property 1 and discussions in the proof. In Lemma 1, the presented density (20) outlines a formula for multivariate bell-curve distributions dependent on the value of q. As q shifts from values approaching 11 from above to values approaching 1+2n1+ 2n form below, the resulting density transitions from Gaussian through a scaled version of the multivariate t−t-distribution to Cauchy and beyond. This density explicitly details all parameters, enabling the application of the maximum likelihood principle and facilitating the use of maximum likelihood estimation in modeling correlated data performed in Section 5. In the context of (20), the characterization matrix 11, denoted by Σ , can undergo modifications to include the degree of freedom parameter m≔2q−1−nm 2q-1-n 58; specifically, p(x,q,Λ) p (x;q, ) =1|πΛ|1/2⋅Γ(m2+n2)Γ(m2)⋅(1+⟨x,Λ−1x⟩)11−q, = 1 |π | 12· ( m2+ n2 ) ( m2 )· (1+ x, ^-1x ) 11-q, (31) where Λ ≔mΣ. m . and m m ≔2q−1−n 2q-1-n Below are a few useful remarks related to the qqGaussian distribution and other bell-curve distributions. Remark 2. To incorporate the location parameter μ, (20) and (31) become p(x,μ,q,Σ) p (x;μ,q, ) =1|πΣ|1/2⋅Γ(1q−1)Γ(1q−1−n2)⋅(2q−1−n)−n2 = 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 )· ( 2q-1-n )^- n2 ⋅(1+(2q−1−n)−1⋅⟨x−μ,Σ−1(x−μ)⟩)11−q; · (1+ ( 2q-1-n )^-1· x-μ, ^-1 (x-μ ) ) 11-q; p(x,μ,q,Λ) p (x;μ,q, ) =1|πΛ|1/2⋅Γ(m2+n2)Γ(m2)⋅(1+⟨x−μ,Λ−1(x−μ)⟩)11−q. = 1 |π | 12· ( m2+ n2 ) ( m2 )· (1+ x-μ, ^-1 (x-μ ) ) 11-q. Remark 3. For random vector X∼qGaussian(q,Σ)X qGaussian (q, ), its q−q-covariance is q[XXT]=(1q−1−n2)1−q2|πΣ|1−q2⋅Γ(q−1−n2)/(Γ(1q−1−n2))qΓ(q−1)/(Γ(1q−1))q⋅Σ.E_q [X^T ]= ( 1q-1- n2 ) 1-q2 |π | 1-q2· ( qq-1- n2 )/ ( ( 1q-1- n2 ) )^q ( qq-1 )/ ( ( 1q-1 ) )^q· . (32) proof We note that ∫ℝnpq(x,q,Σ)x _R^np^q (x;q, )dx =∫ℝn(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2)⋅(1+⟨x,Λ−1x⟩)11−q)qx = _R^n ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 )· (1+ x, ^-1x ) 11-q )^qdx =(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q⋅∫ℝn((1+⟨x,Λ−1x⟩)11−q)qx = ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q· _R^n ( (1+ x, ^-1x ) 11-q )^qdx =(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q⋅∫ℝn(1+⟨x,Λ−1x⟩)q1−qx = ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q· _R^n (1+ x, ^-1x ) q1-qdx =(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q|Λ|1/2⋅∫ℝn(1+⟨x,x⟩)q1−qx = ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q | | 12· _R^n (1+ x,x ) q1-qdx =2n(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q|Λ|1/2⋅∫ℝ>0n(1+⟨x,x⟩)q1−qx =2^n ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q | | 12· _R_>0^n (1+ x,x ) q1-qdx =2n(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q|Λ|1/2⋅∫0∞rn−1(1+r2)q1−q =2^n ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q | | 12· _0^∞r^n-1 (1+r^2 ) q1-q ⋅(∏i=1n−2∫0π2sinn−1−i(θi)dθi⋅∫0π21θ)r · ( _i=1^n-2 _0 π2 ^n-1-i( _i)d _i· _0 π21dθ )dr =π2n−1(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q|Λ|1/2⋅(∏i=1n−212Γ(n−i2)πΓ(n−i+12)) =π 2^n-1 ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q | | 12· ( _i=1^n-2 12 ( n-i2 ) π ( n-i+12 ) ) ⋅∫0∞rn−1(1+r2)q1−qr · _0^∞r^n-1 (1+r^2 ) q1-qdr =2πn2(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q|Λ|1/2⋅(Γ(n2))−1⋅∫0∞(r1−n(1+r2)q−1)−1r =2π n2 ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q | | 12· ( ( n2 ) )^-1· _0^∞ (r^1-n (1+r^2 ) qq-1 )^-1dr =2πn2(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q|Λ|1/2⋅(Γ(n2))−1⋅12B(q−1−n2,n2) =2π n2 ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q | | 12· ( ( n2 ) )^-1· 12B ( qq-1- n2, n2 ) =πn2(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q|Λ|1/2⋅(Γ(n2))−1⋅Γ(q−1−n2)Γ(n2)Γ(q−1) =π n2 ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q | | 12· ( ( n2 ) )^-1· ( qq-1- n2 ) ( n2 ) ( qq-1 ) =(1|πΛ|1/2⋅Γ(1q−1)Γ(1q−1−n2))q|πΛ|1/2⋅Γ(q−1−n2)Γ(q−1) = ( 1 |π | 12· ( 1q-1 ) ( 1q-1- n2 ) )^q |π | 12· ( qq-1- n2 ) ( qq-1 ) =|πΛ|1−q2⋅Γ(q−1−n2)/(Γ(1q−1−n2))qΓ(q−1)/(Γ(1q−1))q. = |π | 1-q2· ( qq-1- n2 )/ ( ( 1q-1- n2 ) )^q ( qq-1 )/ ( ( 1q-1 ) )^q. From (13), we then have the following expression for the q−q-variance-covariance matrix q[XXT] _q [X^T ] =∫ℝnxxT⋅pq(x,q,Σ)x = _R^nx^T· p^q (x;q, )dx =∫ℝnpq(x,q,Σ)x⋅Σ = _R^np^q (x;q, )dx· =|πΛ|1−q2⋅Γ(q−1−n2)/(Γ(1q−1−n2))qΓ(q−1)/(Γ(1q−1))q⋅Σ = |π | 1-q2· ( qq-1- n2 )/ ( ( 1q-1- n2 ) )^q ( qq-1 )/ ( ( 1q-1 ) )^q· =m1−q2|πΣ|1−q2⋅Γ(q−1−n2)/(Γ(1q−1−n2))qΓ(q−1)/(Γ(1q−1))q⋅Σ =m 1-q2 |π | 1-q2· ( qq-1- n2 )/ ( ( 1q-1- n2 ) )^q ( qq-1 )/ ( ( 1q-1 ) )^q· =(1q−1−n2)1−q2|πΣ|1−q2Σ⋅Γ(q−1−n2)/(Γ(1q−1−n2))qΓ(q−1)/(Γ(1q−1))q. = ( 1q-1- n2 ) 1-q2 |π | 1-q2 · ( qq-1- n2 )/ ( ( 1q-1- n2 ) )^q ( qq-1 )/ ( ( 1q-1 ) )^q. (33) Remark 4. Multivariate t distributions with degree of freedom of 2q−1−n 2q-1-n and the scale matrix of Σ are qqGaussian with shape parameter q and scale matrix Σ . Remark 5. The multivariate Cauchy distributions 31 in ℝnR^n with scale matrix 12Σ 12 are qqGaussian with shape parameter q=1+2n+1q=1+ 2n+1 and scale matrix Σ . Remark 6. For a random vector X∼qGaussian(q,Σ)X qGaussian (q, ), the variance-covariance matrix exists if and only if q<1+2n+2q<1+ 2n+2; following the same procedure to derive the variance-covariance matrix for multivariate t distribution yields that, if existing, [XXT]=m−2Σ,E [X^T ]= mm-2 , (34) where m=2q−1−nm= 2q-1-n. The remarks above on qqGaussian distributions reveal their flexibility in incorporating a location parameter, μ, and adapting to multivariate contexts through detailed formulas. Remarkably, these distributions bridge with the class of multivariate bell curve distributions including Gaussian, scaled t, and Cauchy distributions under certain conditions on the shape parameter q. The elaboration of q−q-correlation and the q−q-variable-covariance matrix underscores the capability of these distributions to model and understand the intricacies of correlated data effectively. 4 Proximal Conjugate Gradient Algorithm As delineated in Section 5.3, our optimization scenario is predominantly quadratic in nature; therefore, the conjugate gradient approach has potential for fast convergence and numerical stability. This insight forms the basis for our introduction of a proximal conjugate gradient algorithm framework, tailored to navigate the complexities introduced by the nonconvex penalized qqGaussian likelihood function for sparse statistical learning. To lay the groundwork for this discussion, we begin with a overview of relevant concepts in variational and nonsmooth analysis, presented in Section 4.1. The results presented in Section 4.1 can be found in recent textbooks on variational and nonsmooth analysis, such as 46; 9; 35; 36; 37; 3. 4.1 A Review on Variational and Nonsmooth Analysis Let k,αHC^k, _H with k∈ℕ≥0k _≥ 0 and αH∈[0,1] _H∈ [0,1 ] denote the function space such that ∀F∈k,αH∀ F ^k, _H, F is kkth continuously differentiable, and DkFD^kF is globally Hölder continuous with exponent αH _H; clearly, when αH=1 _H=1, DkFD^kF is globally Lipschitz continuous. In this subsection, we will state the results from variational and nonsmooth analysis related to the following optimization problem: minx∈ℝp+1f(x)≔g(x)+h(x), \ _x ^p+1f (x ) g (x )+h (x ), (35) where f∈0,0(ℝp+1,ℝ)f ^0,0 (R^p+1,R ) is a locally-Lipschitz proper function, g∈1,1(ℝp+1,ℝ)g ^1,1 (R^p+1,R ) is globally L∇g−L_∇ g-smooth and possibly nonconvex, and h∈0,0(ℝp+1,ℝ)h ^0,0 (R^p+1,R ) is a convex locally-Lipschitz function, possibly nonsmooth. The globally Lipschitz property of ∇g∇ g can be alternatively addressed by carrying out the optimization over a compact set. In such scenarios, given that ∇g∇ g is locally Lipschitz, it inherently becomes globally Lipschitz when restricted to a compact set. Results from convex analysis suggest that g,hg,h are Clarke regular; thus, f is Clarke regular. The Clarke’s directional derivative, defined by f∘(x,d) f (x;d ) ≔limy→xsupt↘0f(y+td)−f(y)t _y→ x _t 0 f (y+td )-f (y )t =infδ>0sup‖y−x‖≤δ,0<t<δf(y+td)−f(y)t, = _δ>0 _ \|y-x \|≤δ,0<t<δ f (y+td )-f (y )t, exists for all x∈ℝp+1x ^p+1 since f is Clarke regular. The Clarke subdifferential, denoted by ∂∘ _ , is a set-valued mapping defined by ∂∘f(x)≔ϕ∈ℝp+1|∀d∈ℝp+1,⟨ϕ,d⟩≤f∘(x;d). _ f (x ) \φ ^p+1|∀ d ^p+1,\ φ,d ≤ f (x;d ) \. (36) Since f is a locally Lipschitz function, ∀x∈ℝp+1,∂∘f(x)≠∅∀ x ^p+1,\ _ f (x )≠ . Fundamental convex analysis results show that ∀x∈ℝp+1,∂∘f(x)∀ x ^p+1,\ _ f (x ) is compact, convex, and upper-semicontinuous. ∀x,d∈ℝp+1∀ x,d ^p+1, and we also have f∘(x,d)=maxu∈∂∘f(x)⟨u,d‖d‖⟩.f (x;d )= \ _u∈ _ f (x ) u, d \|d \| . (37) Furthermore, (37) is upper-semicontinuous with respect to x. Simple convex geometry results conclude that (v,−1)|v∈∂∘f(x)=Nepi f(x,f(x)), \ (v,-1 )|v∈ _ f (x ) \=N_epi f (x,f (x ) ), (38) where Nepi f(x,f(x))N_epi f (x,f (x ) ) denotes the normal cone to epi fepi f at the point (x,f(x)) (x,f (x ) ). Since g is smooth, ∂∘g(x)=∇g(x) _ g (x )= \∇ g (x ) \ is a singleton. Then ∂∘f(x)=∂∘g(x)+∂∘h(x) _ f (x )= _ g (x )+ _ h (x ), and f∘(x,d) f (x;d ) =maxu∈∂∘f(x)⟨u,d‖d‖⟩ = \ _u∈ _ f (x ) u, d \|d \| =maxu∈(∇g(x)+∂∘h(x))⟨u,d‖d‖⟩ = \ _u∈ (∇ g (x )+ _ h (x ) ) u, d \|d \| =⟨∇g(x),d‖d‖⟩+maxv∈∂∘h(x)⟨v,d‖d‖⟩ = ∇ g (x ), d \|d \| + \ _v∈ _ h (x ) v, d \|d \| (39) =g∘(x,d)+h∘(x,d). =g (x;d )+h (x;d ). Let Mρt(x)≔(t□(12ρ‖⋅‖2))(x)=infy∈ℝp+1t(y)+12ρ‖y−x‖2M_ρt (x ) (t ( 12ρ \|· \|^2 ) ) (x )= _y ^p+1t (y )+ 12ρ \|y-x \|^2 (40) denote the Moreau envelope operator parameterized by ρ∈ℝ>0ρ _>0 applied on an arbitrary proper, lower semi-continuous, locally Lipschitz function t∈0,0(ℝp+1,ℝ)t ^0,0 (R^p+1,R ), where “□ ” denotes the infimal convolution operator. We have that the Moreau envelope is a smoothing operator, specifically, epi t+epi 12ρ‖⋅‖2⊆epi Mρt,epi t+epi 12ρ \|· \|^2 M_ρt, (41) where “epi” denotes the epigraph. Clearly, Mρt(x)≤t(x)M_ρt (x )≤ t (x ), since (0,0)∈epi 12ρ‖⋅‖2 (0,0 ) 12ρ \|· \|^2 implies that epi t=epi t+(0,0)⊆epi t+epi 12ρ‖⋅‖2⊆epi Mρtepi t=epi t+ (0,0 ) t+epi 12ρ \|· \|^2 M_ρt. When t is convex, (41) takes the equal sign; i.e., the infimal convolution becomes the exact infimal convolution. Consider the affine function A(x)≔⟨a,x⟩+b,A (x ) a,x +b, (42) simple algebra shows that the Moreau envelope applied on A is MρA(x)=⟨a,x⟩+b+ρ2‖a‖2=A(x)+ρ2‖a‖2M_ρA (x )= a,x +b+ ρ2 \|a \|^2=A (x )+ ρ2 \|a \|^2 (43) for some a,b∈ℝp+1a,b ^p+1. Moreover, the following affine addition property is often used in proximal algorithms, mainly due to the fact that the epigraph of an affine function is a half-space that the Moreau envelope applied on: Mρ(t+A)(x) M_ρ (t+A ) (x ) =Mρt(x−ρa)+⟨a,x⟩+b−ρ2‖a‖2 =M_ρt (x-ρ a )+ a,x +b- ρ2 \|a \|^2 (44) Let proxρt(x)≔argMρt(x)=argminy∈ℝp+1t(y)+12ρ‖y−x‖2prox_ρ t (x ) M_ρt (x )= y ^p+1 \ t (y )+ 12ρ \|y-x \|^2 (45) denote the proximal operator, a set-valued mapping; we have proxρt=(I+ρ∂∘t)−1prox_ρ t= (I+ρ _ t )^-1 (46) is the resolvent of the Clarke’s subdifferential operator ρ∂∘tρ _ t. For nonsmooth problems, proximal methods are often used. Fundamental convex analysis results show that: 1. the Moreau envelope Mρt(x)M_ρt (x ) is twice differentiable; thus, its gradient ∇Mρt(x)∇ M_ρt (x ) is well-defined. 2. If t is convex, proxρt(x)prox_ρ t (x ) is a singleton. For the sake of parsimony, with a slight abuse of notation, we use proxρtprox_ρ t to represent a function in this case. It follows that both proxρtprox_ρ t and ∇Mρt∇ M_ρt are firmly non-expansive, and Moreau’s decomposition theorem implies that ∇Mρt(x)=ρ−1(x−proxρt(x)).∇ M_ρt (x )=ρ^-1 (x-prox_ρ t (x ) ). (47) The results from variational and nonsmooth analysis in this subsection have laid the foundation for proving the properties discussed in Section 4.2. 4.2 Proximal Conjugate Gradient Framework Proximal methods are powerful optimization techniques and are particularly adept at handling problems characterized by sparsity, which usually leads to an optimization problem that is nonsmooth 40. Proximal algorithm tends to outperform other methods by far for nonsmooth problems 64; 32. On another ground, Krylov subspace methods represent a cornerstone of numerical analysis, providing a powerful framework for solving large-scale optimization problems efficiently 49. Krylov subspace methods exhibit a remarkable property of convergence acceleration and vastly improved numerical stability, making them indispensable tools in the numerical analyst’s toolkit. Having reviewed the related results from variational and non-smooth analysis in Section 4.1, we are ready to introduce our main optimization framework to combine proximal methods and conjugate gradient together. The essence of proximal algorithms lies upon the Moreau envelope’s smoothing on the objective function. Indeed, proximal methods minimize MρfM_ρf instead of f, thus avoiding nonsmoothness since MρfM_ρf is a smooth function. In this view, proximal algorithms are, in fact, minimizing the Moreau envelope of the objective function. Thus, a wide class of numerical optimization algorithms can easily have their proximal version. Among those, conjugate gradients, a type of Krylov subspace method, are the state-of-the-art methods in smooth optimization due to their computational and memory efficiency, scalability, and numerical stability. Prior to introducing our proximal conjugate gradient update framework, we will first show the equivalency of the optimization problem to minimize (35) and its the Moreau envelope. In nonconvex optimization, the main task for numerical optimization is to find a Clarke stationary point of the objective function, for which we show in Theorem 7 that the set of Clarke stationary point of f is identical to that of MρfM_ρf for ρ∈(0,L∇g−1)ρ∈ (0,L_∇ g^-1 ). Lemma 7. ∀x¯∈ℝp+1,ρ∈(0,L∇g−1)∀ x ^p+1,ρ∈ (0,L_∇ g^-1 ), 0∈∂∘f(x¯)⇔∇Mρf(x¯)=0,0∈ _ f ( x ) ∇ M_ρf ( x )=0, (48) proof Consider arbitrary x∈ℝp+1,ρ∈(0,L∇g−1)x ^p+1,\ ρ∈ (0,L_∇ g^-1 ). As discussed previously, the gradient of the Moreau envelope ∇Mρf(x)=ρ−1(x−proxρf(x))∇ M_ρf (x )=ρ^-1 (x-prox_ρ f (x ) ) implies that proxρf(x)=x−ρ∇Mρf(x),prox_ρ f (x )=x-ρ∇ M_ρf (x ), (49) which implies the following first-order (necessary) optimality condition for Clarke’s stationary point: 0∈ρ−1(x−ρ∇Mρf(x)−x)+∂∘f(x−ρ∇Mρf(x)).0∈ρ^-1 (x-ρ∇ M_ρf (x )-x )+ _ f (x-ρ∇ M_ρf (x ) ). (50) The relation above is simplified to ∇Mρf(x)∈∂∘f(x−ρ∇Mρf(x))=∇g(x−ρ∇Mρf(x))+∂∘h(x−ρ∇Mρf(x)).∇ M_ρf (x )∈ _ f (x-ρ∇ M_ρf (x ) )=∇ g (x-ρ∇ M_ρf (x ) )+ _ h (x-ρ∇ M_ρf (x ) ). (51) Consider arbitrary x¯∈ℝp+1,ρ∈(0,L∇g−1) x ^p+1,ρ∈ (0,L_∇ g^-1 ). “⇒ ” of (48): Let 0∈∂∘f(x¯)=∇g(x¯)+∂∘h(x¯)0∈ _ f ( x )=∇ g ( x )+ _ h ( x ); i.e., x¯ x is a Clarke stationary point of f. Then −∇g(x¯)∈∂∘h(x¯)-∇ g ( x )∈ _ h ( x ). Since h is convex, (51) implies that ⟨−∇g(x¯)−(∇Mρf(x¯)−∇g(x¯−ρ∇Mρf(x¯))),ρ∇Mρf(x¯)⟩≥0. -∇ g ( x )- (∇ M_ρf ( x )-∇ g ( x-ρ∇ M_ρf ( x ) ) ),ρ∇ M_ρf ( x ) ≥ 0. (52) Simplification gives ⟨∇g(x¯−ρ∇Mρf(x¯))−∇g(x¯),∇Mρf(x¯)⟩≥‖∇Mρf(x¯)‖2. ∇ g ( x-ρ∇ M_ρf ( x ) )-∇ g ( x ),∇ M_ρf ( x ) ≥ \|∇ M_ρf ( x ) \|^2. (53) By Cauchy-Schwartz inequality, ⟨∇g(x¯−ρ∇Mρf(x¯))−∇g(x¯),∇Mρf(x¯)⟩≤L∇g⋅ρ‖∇Mρf(x¯)‖2. ∇ g ( x-ρ∇ M_ρf ( x ) )-∇ g ( x ),∇ M_ρf ( x ) ≤ L_∇ g·ρ \|∇ M_ρf ( x ) \|^2. (54) Since ρ<L∇g−1ρ<L_∇ g^-1, (53) and (54) imply that ‖∇Mρf(x¯)‖2≤⟨∇g(x¯−ρ∇Mρf(x¯))−∇g(x¯),∇Mρf(x¯)⟩<‖∇Mρf(x¯)‖2, \|∇ M_ρf ( x ) \|^2≤ ∇ g ( x-ρ∇ M_ρf ( x ) )-∇ g ( x ),∇ M_ρf ( x ) < \|∇ M_ρf ( x ) \|^2, (55) which implies that ∇Mρf(x¯)=0;∇ M_ρf ( x )=0; (56) i.e., x¯ x is the stationary point of MρfM_ρf, hence a Clarke’s stationary point. “⇐ ” of (48): Let ∇fρ(x¯)=0∇ f_ρ ( x )=0; i.e. x¯ x is a stationary point of MρfM_ρf. It follows directly from (51) that 0=∇Mρf(x¯)∈∂∘f(x¯−ρ∇Mρf(x¯))=∂∘f(x¯);0=∇ M_ρf ( x )∈ _ f ( x-ρ∇ M_ρf ( x ) )= _ f ( x ); (57) i.e., x¯ x is a Clarke stationary point of f. The vast majority of optimization algorithms for smooth objective functions require Lipschitz continuity of the objective function. Thus, we are to propose the following Lemma to show the Lipschitz continuity of the gradient of the Moreau envelope of f. Lemma 8. ∀ρ∈(0,L∇g−1)∀ρ∈ (0,L_∇ g^-1 ), ∃L∇Mρf∈ℝ>0∃ L_∇ M_ρf _>0 such that ∀x,y∈ℝp+1,‖∇Mρf(x)−∇Mρf(y)‖≤L∇Mρf‖x−y‖.∀ x,y ^p+1,\ \|∇ M_ρf (x )-∇ M_ρf (y ) \|≤ L_∇ M_ρf \|x-y \|. (58) proof Consider arbitrary x,y∈ℝp+1x,y ^p+1. From (51), since h is convex, ⟨∇Mρf(x)−∇g(x−ρ∇Mρf(x))−(∇Mρf(y)−∇g(y−ρ∇Mρf(y))),x−ρ∇Mρf(x)−(y−ρ∇Mρf(y))⟩≥0. ∇ M_ρf (x )-∇ g (x-ρ∇ M_ρf (x ) )- (∇ M_ρf (y )-∇ g (y-ρ∇ M_ρf (y ) ) ),x-ρ∇ M_ρf (x )- (y-ρ∇ M_ρf (y ) ) ≥ 0. (59) Simplification gives ⟨∇Mρf(x)−∇Mρf(y)−(∇g(x−ρ∇Mρf(x))−∇g(y−ρ∇Mρf(y))),x−y−ρ(∇Mρf(x)−∇Mρf(y))⟩≥0. ∇ M_ρf (x )-∇ M_ρf (y )- (∇ g (x-ρ∇ M_ρf (x ) )-∇ g (y-ρ∇ M_ρf (y ) ) ),x-y-ρ (∇ M_ρf (x )-∇ M_ρf (y ) ) ≥ 0. (60) Let δ∇Mρf≔∇Mρf(x)−∇Mρf(y) _∇ M_ρf ∇ M_ρf (x )-∇ M_ρf (y ), δ∇g≔∇g(x−ρ∇fρ(x))−∇g(y−ρ∇fρ(y)) _∇ g ∇ g (x-ρ∇ f_ρ (x ) )-∇ g (y-ρ∇ f_ρ (y ) ), and δx,y≔x−y _x,y x-y, then 0 0 ≤⟨δ∇Mρf−δ∇g,δx,y−ρδ∇Mρf⟩ ≤ _∇ M_ρf- _∇ g, _x,y-ρ _∇ M_ρf =−ρ‖δ∇Mρf‖2+ρ⟨δ∇g,δ∇Mρf⟩+⟨δ∇Mρf,δx,y⟩−⟨δ∇g,δx,y⟩ =-ρ \| _∇ M_ρf \|^2+ρ _∇ g, _∇ M_ρf + _∇ M_ρf, _x,y - _∇ g, _x,y ≤−ρ‖δ∇Mρf‖2+ρ‖δ∇g‖⋅‖δ∇Mρf‖+‖δ∇Mρf‖⋅‖δx,y‖+‖δ∇g‖⋅‖δx,y‖ ≤-ρ \| _∇ M_ρf \|^2+ρ \| _∇ g \|· \| _∇ M_ρf \|+ \| _∇ M_ρf \|· \| _x,y \|+ \| _∇ g \|· \| _x,y \| ≤−ρ‖δ∇Mρf‖2+ρL∇g(‖δx,y‖+ρ‖δ∇Mρf‖)⋅‖δ∇Mρf‖ ≤-ρ \| _∇ M_ρf \|^2+ρ L_∇ g ( \| _x,y \|+ρ \| _∇ M_ρf \| )· \| _∇ M_ρf \| +‖δ∇Mρf‖⋅‖δx,y‖+L∇g(‖δx,y‖+ρ‖δ∇Mρf‖)⋅‖δx,y‖ + \| _∇ M_ρf \|· \| _x,y \|+L_∇ g ( \| _x,y \|+ρ \| _∇ M_ρf \| )· \| _x,y \| Simplification of the above inequality gives ‖δ∇Mρf‖≤2Lgρ+1+8Lgρ+12ρ(1−L∇gρ)‖δx,y‖; \| _∇ M_ρf \|≤ 2L_gρ+1+ 8L_gρ+12ρ (1-L_∇ gρ ) \| _x,y \|; (61) i.e., ‖∇Mρf(x)−∇Mρf(y)‖≤L∇Mρf‖x−y‖, \|∇ M_ρf (x )-∇ M_ρf (y ) \|≤ L_∇ M_ρf \|x-y \|, (62) where L∇Mρf≔2Lgρ+1+8Lgρ+12ρ(1−L∇gρ)>0.L_∇ M_ρf 2L_gρ+1+ 8L_gρ+12ρ (1-L_∇ gρ )>0. (63) Following this idea, we introduce our proximal conjugate gradient framework in Algorithm 1. 1: Input: A fixed value of ρ∈(0,ρ−1)ρ∈ (0,ρ^-1 ) 2: Calculate the gradient of the Moreau envelope: s(k)≔∇Mρf(x(k))s (k ) ∇ M_ρf (x (k ) ) 3: d(k)≔−s(k)+β(k)⋅d(k−1)d (k ) -s (k )+β (k )· d (k-1 ) 4: Line search to find α(k)α (k ) for the update x(k+1)≔x(k)+α(k)d(k)x (k+1 ) x (k )+α (k )d (k ) 5: Update x(k+1)≔x(k)+α(k)d(k)x (k+1 ) x (k )+α (k )d (k ) Algorithm 1 Proximal Point Algorithm In the above algorithm, β(k)β (k ) is the conjugate parameter. The significant meaning of Algorithm 1 is that for any global convergent numerical method to find the equilibria of a globally Lipschitz flow, which generally include the global convergent first-order methods, Algorithm 1 can transform such a method to a proximal counterpart. For some objective functions, the gradient of the Moreau envelope can be calculated directly. However, calculation for the Moreau envelope’s gradient is not tractable for many objective functions whose smooth component g is of complicated form. Motivated by this, we further consider the following the Moreau envelope of the objective function with linearized g, such linearization step is frequently used in proximal algorithms for statistical sparse learning problems (e.g., 38; 19; 63). Consider the linearized surrogate of (35), the locally Lipschitz function f~∈0,0(ℝp+1,ℝ) f ^0,0 (R^p+1,R ), defined by f~(x,u) f (x;u ) ≔⟨u,x⟩+h(x) u,x +h (x ) (64) proxρf~(x,u) _ρ f (x;u ) =argminy∈ℝp+1⟨u,y⟩+12ρ‖y−x‖2+h(y) = y ^p+1 \ \ u,y + 12ρ \|y-x \|^2+h (y ) \ (65) ∇xMρf~(x,u) _xM_ρ f (x;u ) =ρ−1(x−proxρf~(x,u)) =ρ^-1 (x-prox_ρ f (x;u ) ) (66) proxρf~(x,u)prox_ρ f (x;u ) is the proximal operator applied on f~ f, and ∇xMρf~(x,u) _xM_ρ f (x;u ) is the gradient of the Moreau envelope of f~ f. The linearization term ⟨u,x⟩ u,x in (64) depends on u. Recognize that f~(x,u) f (x;u ) is linearizing the nonconvex smooth component g in (101) when u=∇g(x)u=∇ g (x ). We establish several definitions for subsequent utilization. Define the mapping g~ρ=I−ρ∇g∈0,0(ℝp+1,ℝp+1) g_ρ=I-ρ∇ g ^0,0 (R^p+1,R^p+1 ) for some ρ∈(0,L∇g−1)ρ∈ (0,L_∇ g^-1 ), the locally Lipschitz property of g~ρ g_ρ follows from g∈1,1g ^1,1; i.e., g~ρ(x)≔x−ρ∇g(x). g_ρ (x ) x-ρ∇ g (x ). The following Lemma identifies some fundamental property of g~ρ g_ρ. Lemma 9. g~ρ g_ρ is a bijective from ℝp+1R^p+1 to ℝp+1R^p+1, and g~ρ−1 g_ρ^-1 is globally Lipschitz with constant (1−ρL∇g)−1 (1-ρ L_∇ g )^-1. proof Injectivity proof: Consider arbitrary x1,x2∈ℝp+1x_1,x_2 ^p+1. Since ρ∈(0,L∇g−1)ρ∈ (0,L_∇ g^-1 ), x1−ρ∇g(x1)=x2−ρ∇g(x2)x_1-ρ∇ g (x_1 )=x_2-ρ∇ g (x_2 ) implies that ‖x1−x2‖=ρ‖∇g(x1)−∇g(x2)‖≤ρL∇g‖x1−x2‖<‖x1−x2‖, \|x_1-x_2 \|=ρ \|∇ g (x_1 )-∇ g (x_2 ) \|≤ρ L_∇ g \|x_1-x_2 \|< \|x_1-x_2 \|, (67) hence x1=x2x_1=x_2. This shows that g~ρ g_ρ is a injective mapping. Surjectivity proof: Consider arbitrary y1,y2∈ℝp+1y_1,y_2 ^p+1. Consider arbitrary z∈ℝp+1z ^p+1. Define mapping (y)≔z+ρ∇g(y)T (y ) z+ρ∇ g (y ), then ‖(y1)−(y2)‖ \|T (y_1 )-T (y_2 ) \| =‖z+ρ∇g(y1)−(z+ρ∇g(y2))‖ = \|z+ρ∇ g (y_1 )- (z+ρ∇ g (y_2 ) ) \| =ρ‖∇g(y1)−∇g(y2)‖ =ρ \|∇ g (y_1 )-∇ g (y_2 ) \| ≤ρL∇g‖y1−y2‖ ≤ρ L_∇ g \|y_1-y_2 \| <‖y1−y2‖. < \|y_1-y_2 \|. Thus, T is a contraction mapping, since ℝp+1R^p+1 equipped with Euclidean topology is a Banach space, by Banach fixed point theorem, T has a fixed point; i.e., ∃y∈ℝp+1∃ y ^p+1 such that y=z+ρ∇g(y)y=z+ρ∇ g (y ), or equivalently, g~ρ(y)=y−ρ∇g(y)=z g_ρ (y )=y-ρ∇ g (y )=z. Thus, ℝp+1⊆g~ρ(ℝp+1)R^p+1 g_ρ (R^p+1 ). Globally Lipschitz constant derivation for inverse map: Since ∇g∇ g is globally L∇g−L_∇ g-Lipschitz, ‖g~ρ(y1)−g~ρ(y2)‖ \| g_ρ (y_1 )- g_ρ (y_2 ) \| =‖y1−ρ∇g(y1)−(y2−ρ∇g(y2))‖ = \|y_1-ρ∇ g (y_1 )- (y_2-ρ∇ g (y_2 ) ) \| =‖y1−y2−ρ(∇g(y1)−∇g(y2))‖ = \|y_1-y_2-ρ (∇ g (y_1 )-∇ g (y_2 ) ) \| ≥|‖y1−y2‖−‖ρ(∇g(y1)−∇g(y2))‖| ≥ | \|y_1-y_2 \|- \|ρ (∇ g (y_1 )-∇ g (y_2 ) ) \| | =|‖y1−y2‖−ρ‖∇g(y1)−∇g(y2)‖| = | \|y_1-y_2 \|-ρ \|∇ g (y_1 )-∇ g (y_2 ) \| | =‖y1−y2‖−ρ‖∇g(y1)−∇g(y2)‖ = \|y_1-y_2 \|-ρ \|∇ g (y_1 )-∇ g (y_2 ) \| (68) ≥(1−ρL∇g)‖y1−y2‖ ≥ (1-ρ L_∇ g ) \|y_1-y_2 \| where (68) is due to the fact that ρ‖∇g(y1)−∇g(y2)‖≤ρL∇g‖y1−y2‖<‖y1−y2‖.ρ \|∇ g (y_1 )-∇ g (y_2 ) \|≤ρ L_∇ g \|y_1-y_2 \|< \|y_1-y_2 \|. (69) Since g~ρ g_ρ is surjective, consider arbitrary z1,z2∈ℝp+1z_1,z_2 ^p+1 let y1≔g~ρ−1(z1)y_1 g_ρ^-1 (z_1 ) and y2≔g~ρ−1(z2)y_2 g_ρ^-1 (z_2 ), then ‖g~ρ−1(z1)−g~ρ−1(z2)‖≤(1−ρL∇g)−1‖z1−z2‖. \| g_ρ^-1 (z_1 )- g_ρ^-1 (z_2 ) \|≤ (1-ρ L_∇ g )^-1 \|z_1-z_2 \|. (70) Define ρf~(x)≔∇xMρf~(x,u)G_ρ f (x ) _xM_ρ f (x;u ) (71) with u=∇g(x)u=∇ g (x ); i.e., ρf~(x)G_ρ f (x ) is the gradient of the Moreau envelope of f~ f. Similarly to Lemma 7 and 8, we are to prove that the set of Clarke’s stationary of (101) is identical to the set x¯∈ℝp+1|ρf~(x¯)=0 \ x ^p+1|G_ρ f ( x )=0 \ in Lemma 10, and then we are to show that (71) is globally Lipschitz in Lemma 11. Lemma 10. ∀x¯∈ℝp+1,ρ∈ℝ>0∀ x ^p+1,ρ _>0, 0∈∂∘f(x¯)⇔ρf~(x¯)=0.0∈ _ f ( x ) _ρ f ( x )=0. (72) proof Consider arbitrary x∈ℝp+1x ^p+1. The f~ f is convex since it is a sum of convex function h and a linear mapping of x, which is convex. ρf~(x) _ρ f (x ) =ρ−1(x−proxρf~(x,∇g(x))) =ρ^-1 (x-prox_ρ f (x;∇ g (x ) ) ) (73) =ρ−1(x−proxρh(x−ρ∇g(x))) =ρ^-1 (x-prox_ρ h (x-ρ∇ g (x ) ) ) (74) =∇g(x)+ρ−1(x−ρ∇g(x)−proxρh(x−ρ∇g(x))) =∇ g (x )+ρ^-1 (x-ρ∇ g (x )-prox_ρ h (x-ρ∇ g (x ) ) ) =∇g(x)+(∇Mρh)∘g~ρ(x) =∇ g (x )+ (∇ M_ρh ) g_ρ (x ) (75) (74) is due to the affine addition property of proximal mapping. From (64) and (73), ρf~(x) _ρ f (x ) =ρ−1(x−proxρf~(x,∇g(x))) =ρ^-1 (x-prox_ρ f (x,∇ g (x ) ) ) ⟹proxρf~(x,∇g(x)) _ρ f (x,∇ g (x ) ) =x−ρ⋅ρf~(x) =x-ρ·G_ρ f (x ) ⟹0 0 ∈ρ−1(x−ρ⋅ρf~(x)−x)+∂∘f~(x−ρ⋅ρf~(x)) ∈ρ^-1 (x-ρ·G_ρ f (x )-x )+ _ f (x-ρ·G_ρ f (x ) ) ⟹0 0 ∈−ρf~(x)+∇g(x)+∂∘h(x−ρ⋅ρf~(x)) ∈-G_ρ f (x )+∇ g (x )+ _ h (x-ρ·G_ρ f (x ) ) ⟹ρf~(x) _ρ f (x ) ∈∇g(x)+∂∘h(x−ρ⋅ρf~(x)) ∈∇ g (x )+ _ h (x-ρ·G_ρ f (x ) ) (76) ⟹ρf~(x)−∇g(x) _ρ f (x )-∇ g (x ) ∈∂∘h(x−ρ⋅ρf~(x)). ∈ _ h (x-ρ·G_ρ f (x ) ). Thus, since h is convex, ∀v∈∂∘h(x)∀ v∈ _ h (x ), ⟨ρf~(x)−∇g(x)−v,x−ρ⋅ρf~(x)−x⟩ _ρ f (x )-∇ g (x )-v,x-ρ·G_ρ f (x )-x ≥0 ≥ 0 ⟹⟨ρf~(x)−∇g(x)−v,ρf~(x)⟩ _ρ f (x )-∇ g (x )-v,G_ρ f (x ) ≤0 ≤ 0 ⟹‖ρf~(x)‖2 \|G_ρ f (x ) \|^2 ≤⟨∇g(x)+v,ρf~(x)⟩ ≤ ∇ g (x )+v,G_ρ f (x ) (77) ≤‖∇g(x)+v‖⋅‖ρf~(x)‖ ≤ \|∇ g (x )+v \|· \|G_ρ f (x ) \| ⟹‖ρf~(x)‖ \|G_ρ f (x ) \| ≤‖∇g(x)+v‖, ≤ \|∇ g (x )+v \|, (78) provided that ‖ρf~(x)‖≠0 \|G_ρ f (x ) \|≠ 0. Basic results on the Moreau envelope shows that ρf~(x)=0G_ρ f (x )=0 implies that x is a Clarke stationary point of f~(x) f (x ). Now we are proceed to prove (72): “⇒ ”: Consider arbitrary x¯∈ℝp+1 x ^p+1 and ρ∈ℝ>0ρ _>0. Let 0∈∂∘f(x¯)=∇g(x¯)+∂∘h(x¯)0∈ _ f ( x )=∇ g ( x )+ _ h ( x ); i.e., x¯ x is a Clarke stationary point of f. Then ∃v∈∂∘h(x¯)∃ v∈ _ h ( x ) such that ∇g(x¯)+v=0∇ g ( x )+v=0. (78) implies that ‖ρf~(x¯)‖≤‖∇g(x¯)+v‖=0. \|G_ρ f ( x ) \|≤ \|∇ g ( x )+v \|=0. (79) Thus, ρf~(x¯)=0G_ρ f ( x )=0. “⇐ ”: Consider arbitrary x¯∈ℝp+1 x ^p+1. Let ρf~(x¯)=0G_ρ f ( x )=0; i.e., ρf~(x¯)=0G_ρ f ( x )=0 is stationary. (76) implies that 0=ρf~(x¯)∈∇g(x¯)+∂∘h(x¯−ρ⋅ρf~(x¯))=∇g(x¯)+∂∘h(x¯)=∂∘f(x¯).0=G_ρ f ( x )∈∇ g ( x )+ _ h ( x-ρ·G_ρ f ( x ) )=∇ g ( x )+ _ h ( x )= _ f ( x ). (80) Thus, x¯ x is a Clarke stationary point of f. Lemma 11. ∀ρ∈ℝ>0∀ρ _>0, ∃Lρf~∈ℝ>0∃ L_G_ρ f _>0 such that ∀x,y∈ℝp+1,‖ρf~(x)−ρf~(y)‖≤Lρf~‖x−y‖.∀ x,y ^p+1,\ \|G_ρ f (x )-G_ρ f (y ) \|≤ L_G_ρ f \|x-y \|. (81) proof Consider arbitrary x,y∈ℝp+1x,y ^p+1 and ρ∈ℝ>0ρ _>0. Let u≔∇g(x)u ∇ g (x ) and v≔∇g(y)v ∇ g (y ), ‖ρf~(x)−ρf~(y)‖ \|G_ρ f (x )-G_ρ f (y ) \| =‖∇xMρf~(x,u)−∇yMρf~(y,v)‖ = \| _xM_ρ f (x;u )- _yM_ρ f (y;v ) \| =‖∇xMρf~(x,u)−∇xMρf~(x,v)+∇xMρf~(x,v)−∇yMρf~(y,v)‖ = \| _xM_ρ f (x;u )- _xM_ρ f (x;v )+ _xM_ρ f (x;v )- _yM_ρ f (y;v ) \| ≤‖∇xMρf~(x,u)−∇xMρf~(x,v)‖+‖∇xMρf~(x,v)−∇yMρf~(y,v)‖ ≤ \| _xM_ρ f (x;u )- _xM_ρ f (x;v ) \|+ \| _xM_ρ f (x;v )- _yM_ρ f (y;v ) \| ≤‖u−v‖+‖∇xMρf~(x,v)−∇yMρf~(y,v)‖ ≤ \|u-v \|+ \| _xM_ρ f (x;v )- _yM_ρ f (y;v ) \| (82) ≤‖u−v‖+ρ−1‖x−y‖ ≤ \|u-v \|+ρ^-1 \|x-y \| (83) ≤L∇g‖x−y‖+ρ−1‖x−y‖ ≤ L_∇ g \|x-y \|+ρ^-1 \|x-y \| =(L∇g+ρ−1)‖x−y‖ = (L_∇ g+ρ^-1 ) \|x-y \| (82) is due to Lemma 4 in 19, and (83) is due to f~(⋅,v) f (·;v ) is convex and the fact that the gradient of a convex function’s Moreau envelope is ρ−1−ρ^-1-Lipschitz. Therefore, let Lρf~≔L∇g+ρ−1L_G_ρ f L_∇ g+ρ^-1 (84) and we have ‖ρf~(x)−ρf~(y)‖≤Lρf~‖x−y‖. \|G_ρ f (x )-G_ρ f (y ) \|≤ L_G_ρ f \|x-y \|. (85) Furthermore, (75) suggests that ρf~=∇g+(∇Mρh)∘g~ρ=∇g+(∇Mρh)∘(Id−ρ∇g).G_ρ f=∇ g+ (∇ M_ρh ) g_ρ=∇ g+ (∇ M_ρh ) (Id-ρ∇ g ). (86) Hence, Id−ρρf~ Id- _ρ f =Id−ρ∇g−ρ(∇Mρh)∘(Id−ρ∇g) =Id-ρ∇ g-ρ (∇ M_ρh ) (Id-ρ∇ g ) =g~ρ−ρ(∇Mρh)∘g~ρ = g_ρ-ρ (∇ M_ρh ) g_ρ =(Id−ρ(∇Mρh))∘g~ρ = (Id-ρ (∇ M_ρh ) ) g_ρ (87) =g~ρ−1∘(g~ρ∘(Id−ρ(∇Mρh)))∘g~ρ = g_ρ^-1 ( g_ρ (Id-ρ (∇ M_ρh ) ) ) g_ρ (88) shows that g~ρ−1∘(g~ρ∘(Id−ρ(∇Mρh)))∘g~ρ:ℝp+1↦ℝp+1 g_ρ^-1 ( g_ρ (Id-ρ (∇ M_ρh ) ) ) g_ρ:R^p+1 ^p+1 equals to Id−ρρf~:ℝp+1↦ℝp+1Id- _ρ f:R^p+1 ^p+1. Since g~ρ g_ρ ibijectivee, and that g~ρ g_ρ and g~ρ−1 g_ρ^-1 are continuous due to the globally Lipschitz property from Lemma 9, g~ρ g_ρ is a homeomorphism. Hence, Id−ρρf~Id- _ρ f and g~ρ∘(Id−ρ(∇Mρh)) g_ρ (Id-ρ (∇ M_ρh ) ) are topologically equivalent mappings via the homeomorphism g~ρ g_ρ. Lemma 11 implies that ρf~G_ρ f is globally Lipschitz, which sufficiently implies by the Cauchy-Lipschitz theorem that the differential equation x˙≔dxdt=ρf~(x) x dxdt=G_ρ f (x ) (89) has a unique solution for any given initial value condition. Thus, ρf~G_ρ f generates a unique flow under a given initial value condition. The operator equations presented above can be understood as demonstrating how ρf~G_ρ f functions analogously to a gradient operator. Specifically, Id−ρρf~Id- _ρ f represents executing a descent operation in the −ρf~-G_ρ f direction with a step size ρ. Similarly, Id−ρ(∇Mρh)Id-ρ (∇ M_ρh ) represents a single gradient descent step with ρ as the step size with objective function MρhM_ρh, the Moreau envelope of h; while g~ρ=Id−ρ∇g g_ρ=Id-ρ∇ g reflects a gradient descent step with objective function g, again with ρ as the step size. Equation (87) elucidates that a descent in the −ρf~-G_ρ f direction is identical to first performing a one-step gradient descent on g, followed by MρhM_ρh; or performing gradient descents in a converse order yields a topologically equivalence via the homeomorphism g~ρ g_ρ, as shown in (88). In short summary, the approach based on linearization of the smooth term and the Moreau envelope enables us to build equivalence between identifying Clarke stationary points of the original nonsmooth objective function (101) and finding equilibria of the (unique) flow generated by ρf~G_ρ f, as demonstrated in Lemma 10. The task of finding equilibria within a globally Lipschitz continuous flow, such as the ρf~G_ρ f flow, is well explored within mathematics, particularly in the realms of dynamical systems and numerical analysis (see, for example, 44; 2; 33; 27; 24). Cauchy-Lipschitz theorem establishes the uniqueness of solutions to initial value problems for globally Lipschitz continuous flows; while the existence of equilibria is a direct result of Brouwer fixed-point theorem. Numerical methods for dynamical systems, including methods for finding equilibria of the flow, are largely based on this uniqueness result. This is one reason that the vast majority of numerical methods in the context of dynamical systems require the flow to be globally Lipschitz. It is important to note that these numerical strategies, widely applied across dynamical systems, do not hinge on the flow being derived from a conservative field. As such, the process of formulating a potential function for ρf~G_ρ f is not a prerequisite for employing numerical techniques to determine its equilibria. This perspective underscores the versatility of numerical methods in dynamical systems in finding the equilibria of flows, regardless of the explicit existence of a potential function, a stance corroborated by various sources in the literature 44; 2; 33; 27; 24; 47; 45. In this view, the construction of a potential function for ρf~G_ρ f is generally not necessary when deploying numerical analysis methods to find its equilibria. In the context of nonlinear conjugate gradient algorithms for optimization, achieving global convergence on nonconvex objective functions that are globally Lipschitz-smooth implies that such methods can reliably find equilibria within the corresponding flow dynamics 47; 45. These algorithms typically incorporate a line search step, which may use a surrogate objective function instead of the original. This surrogate can be a constructed potential, Lyapunov, or energy function, offering flexibility when finding the potential function for ρf~G_ρ f poses challenges 47; 10; 54. When it is feasible to construct a potential function whose gradient with respect to x is ρf~G_ρ f, the associated objective function and its gradient become more manageable, allowing for direct global convergence arguments. If constructing a potential function with respect to x for the (∇Mρh)∘g~ρ(x) (∇ M_ρh ) g_ρ (x ) term in (75) or ∇g∘(Id−ρ(∇Mρh))∇ g (Id-ρ (∇ M_ρh ) ) in (88) is tractable, the objective function with gradient being (75) or g~ρ∘(Id−ρ(∇Mρh)) g_ρ (Id-ρ (∇ M_ρh ) ) can hence be easily constructed. Thus, arguments for global convergence for methods based on the objective function and its gradient directly follow to prove the global convergence of the numerical optimization algorithm when applied to the constructed potential function for ρf~G_ρ f. We remark that Id−ρρf~Id- _ρ f and g~ρ∘(Id−ρ(∇Mρh)) g_ρ (Id-ρ (∇ M_ρh ) ) generate two topologically equivalent flows via homeomorphism g~ρ g_ρ; thus, their equilibria can be transformed by g~ρ g_ρ and share the same stability. In the context of numerical optimization, this implies that a fixed point x¯ x for the mapping Id−ρρf~Id- _ρ f corresponds bijectively to a fixed point g~ρ(x¯) g_ρ ( x ) for g~ρ∘(Id−ρ(∇Mρh)) g_ρ (Id-ρ (∇ M_ρh ) ). Characterized by the first-order optimality condition in optimization of smooth functions, or equivalently, the stationary condition in dynamical system, ρf~(x¯)=∇g(x¯)+(∇Mρh)∘g~ρ(x¯)=0⇔∇Mρh(g~ρ(x¯))+∇g(g~ρ(x¯)−ρ∇Mρh(g~ρ(x¯)))=0.G_ρ f ( x )=∇ g ( x )+ (∇ M_ρh ) g_ρ ( x )=0 ∇ M_ρh ( g_ρ ( x ) )+∇ g ( g_ρ ( x )-ρ∇ M_ρh ( g_ρ ( x ) ) )=0. (90) This approach is practical because the literature on first-order numerical optimization techniques frequently includes proofs of global convergence for methods that depend on the objective function and its gradient (for example, see 17; 43; 25; 13; 22). Alternatively, construction of a potential function for ρf~G_ρ f is often not necessary due to the fact that fixed-point methods finding equilibria for a flow mostly establish convergence properties based on Banach fixed point theorem. This theorem guarantees convergence through intrinsic flow properties, obviating the need for a potential function Burden2001; 2; 1. Conventionally, the use of line search based on the objective function and its gradient has been applied in some numerical methods to ensure global convergence. However, with the rapid growth of research in high–dimensional statistical machine learning and large-scale optimization, evaluations of the objective function often proven to be inefficient. Consequently, recent years have seen the exploration of two main alternatives. For instance, two different types of approaches for global convergent nonlinear conjugate gradient methods have been proposed without the conventional objective function-based line search procedure. One type of approach ensures global convergence by utilizing a line search mechanism that depends only on the nonlinear equation that generates the flow 16; 52; 53; 28; that is, the gradient function for smooth optimization, or ρf~G_ρ f in our case. As an example, under the smoothness assumption, the first-order optimality condition for an exact line search often solves for α with the current value x(k)x (k ) and the search direction d(k)d (k ) from ⟨ρf~(x(k)+α⋅d(k)),d(k)⟩=0 _ρ f (x (k )+α· d (k ) ),d (k ) =0, an equation dependent only on ρf~G_ρ f but not any surrogate objective function. From a practical perspective, this one-dimensional root finding problem can be carried out efficiently using the Brent root finding algorithm 7. The other approach suggests achieving global convergence either without the need for line search 51; 8; 55; 62; 61; 65 or by meeting a condition related to the Zoutendijk condition to replace the Wolfe-Powell conditions of sufficient descent (Armijo) and curvature 39. Additionally, in scenarios where the fulfillment of a sufficient descent (Armijo) condition is imperative, the formulation of a surrogate objective function becomes essential. Considering (75), where a surrogate objective is required for the line search phase, it could be formulated as: g(x)+(Mρh)∘g~ρ(x) g (x )+ (M_ρh ) g_ρ (x ) =g(x)+(Mρh)∘g~ρ(x)+constant =g (x )+ (M_ρh ) g_ρ (x )+constant =g(x)+⟨∇g(x),proxρh(x−ρ∇g(x))−x⟩+12ρ‖proxρh(x−ρ∇g(x))−x‖2 =g (x )+ ∇ g (x ),prox_ρ h (x-ρ∇ g (x ) )-x + 12ρ \|prox_ρ h (x-ρ∇ g (x ) )-x \|^2 (91) +h(proxρh(x−ρ∇g(x)))+constant +h (prox_ρ h (x-ρ∇ g (x ) ) )+constant This formulation, denoted as (91), represents a quadratic approximation of g plus the nonsmooth term h, evaluated at proxρh(x−ρ∇g(x))prox_ρ h (x-ρ∇ g (x ) ). This type of formulation has often been used for the line search step in previous studies 4; 29. The addition of the term −⟨∇g(x),x⟩- ∇ g (x ),x acts as a constant in (64), analogous to fixing the value of u as ∇g(x)∇ g (x ) for linearization. This constant term, −⟨∇g(x),x⟩- ∇ g (x ),x , doesn’t alter the gradient of the Moreau envelope (66) or the proximal point (65), serving to frame the quadratic approximation of g(proxρh(x−ρ∇g(x)))g (prox_ρ h (x-ρ∇ g (x ) ) ). Evaluation of proxρhprox_ρ h in (91) is tractable and efficient for many functions, such as the ℓ1 _1 norm commonly encountered in sparse statistical learning can be efficiently computed via the soft-thresholding function. Given that line search rules such as the Wolfe-Powell or Armijo-Goldstein conditions require only the difference in the value of the objective function at two points to decide on the step size, the constant term in (91) can be disregarded. Subsequent global convergence arguments stem from the fixed-point theory analysis of the numerical methods deployed to find the equilibria of the ρf~G_ρ f flow. Another possible surrogate objective function inspired by the quadratic Lyapunov function for the ρf~G_ρ f flow could be 12‖ρf~‖2 12 \|G_ρ f \|^2, attains its minimal value 00 exactly at the ρf~G_ρ f flow’s equilibria. This quadratic approach simplifies evaluation, but it may not offer insights into the potential function’s landscape, potentially limiting the numerical algorithm’s acceleration capabilities if such an algorithm uses the landscape information to ensure the sufficient descent (Armijo) condition. Therefore, formulating the surrogate objective function preserving the landscape of the original objective function as outlined in (91) is preferable. Building on the above discussion, we introduce our practical proximal conjugate gradient framework in Algorithm 2. 1: Input: A fixed value of ρ∈(0,ρ−1)ρ∈ (0,ρ^-1 ) 2: Calculate the proximal value p(k)≔proxρ−1h(x(k)−ρ−1⋅∇g(x(k)))p (k ) _ρ^-1h (x (k )-ρ^-1·∇ g (x (k ) ) ) 3: Calculate ρf~(x(k))G_ρ f (x (k ) ): s(k)≔ρ(x(k)−p(k))s (k ) ρ (x (k )-p (k ) ) 4: d(k)≔−s(k)+β(k)⋅d(k−1)d (k ) -s (k )+β (k )· d (k-1 ) 5: Line search to find α(k)α (k ) for the update x(k+1)≔x(k)+α(k)d(k)x (k+1 ) x (k )+α (k )d (k ), if needed. 6: Update x(k+1)≔x(k)+α(k)d(k)x (k+1 ) x (k )+α (k )d (k ) Algorithm 2 Computationally Tractable Proximal Conjugate Gradient Update Scheme In Algorithm 2, β(k)β (k ) functions as the conjugate parameter. Unlike Algorithm 1, Algorithm 2 facilitates the update process without the need to compute ∇Mρf(x(k))∇ M_ρf (x (k ) ). This adaptation is significantly valuable in practical scenarios, especially in statistical sparse learning challenges characterized by a complicated smooth component g alongside a simple nonsmooth convex component h. In such cases, computing proxρhprox_ρ h is markedly more tractable and efficient than proxρfprox_ρ f. This approach is particularly beneficial for sparse statistical learning issues, where sparsity is commonly induced by an ℓ1 _1 penalty term. 4.3 Proximal Hager-Zhang 22 Conjugate Gradient The nonlinear conjugate gradient method represents the pinnacle of first-order techniques for addressing smooth optimization challenges. Various versions of nonlinear conjugate gradient methods have been introduced, including the Fletcher-Reeves (FR) method 17, the modified Polak-Ribiere-Polyak (PRP+) method 43; 21, the Hestenes-Stiefel (HS) method 25, the Dai-Yuan (DY) method 13, and the Hager-Zhang (HZ) method 22. These versions have all demonstrated global convergence with nonconvex globally Lipschitz-smooth objective functions. Among these, the Hager-Zhang conjugate gradient method is notable for delivering the best numerical performance on large-scale datasets, as indicated in previous research 23. Building on this, having introduced our practical proximal conjugate gradient update mechanism in Algorithm 2, we aim to extend this approach by adapting the smooth Hager-Zhang nonlinear conjugate gradient method to its proximal version in Algorithm 3. 1: Input: Initial point x(0)x (0 ); g∈1,1(ℝp+1,ℝ)g ^1,1 (R^p+1,R ); locally-Lipschitz, convex h∈0,0(ℝp+1,ℝ)h ^0,0 (R^p+1,R ); the smoothing parameter for the Moreau envelope ρ∈(0,ρ−1)ρ∈ (0,ρ^-1 ); k≔0k 0 2: Output: p 3: k+=1k+=1 4: Calculate the gradient for g: g(0)≔∇g(x(0))g (0 ) ∇ g (x (0 ) ) 5: Calculate the proximal value p(0)≔proxρ,h(x(0)−ρ⋅g(0))p (0 ) _ρ,h (x (0 )-ρ· g (0 ) ) 6: Calculate the gradient analog: s(0)≔x(0)−p(0)s (0 ) x (0 )-p (0 ) 7: d(0)≔−s(0)d (0 ) -s (0 ) 8: Perform the line search with d(0)d (0 ) with step size α(0)α (0 ) 9: Update x1≔x(0)+α(0)d(0)x_1 x (0 )+α (0 )d (0 ) 10: while not converged do 11: k+=1k+=1 12: Calculate the gradient for g: g(k)≔∇g(x(k))g (k ) ∇ g (x (k ) ) 13: Calculate the proximal value p(k)≔proxρ,h(x(k)−ρ⋅g(k))p (k ) _ρ,h (x (k )-ρ· g (k ) ) 14: Calculate the gradient analog: s(k)≔x(k)−p(k)s (k ) x (k )-p (k ) 15: d(k)≔−s(k)+β¯(k)⋅d(k−1)d (k ) -s (k )+ β (k )· d (k-1 ) 16: Perform the line search with d(k)d (k ) with step size α(k)α (k ) based on Wolfe-Powell conditions 17: Update x(k+1)≔x(k)+α(k)d(k)x (k+1 ) x (k )+α (k )d (k ) 18: Check for convergence 19: return p(k)p (k ) Algorithm 3 Proximal Hager-Zhang 22 Conjugate Gradient In Algorithm 3, Hager-Zhang’s conjugate parameter β¯(k) β (k ) is defined as 22: y(k) y (k ) ≔s(k+1)−s(k) s (k+1 )-s (k ) β(k) β (k ) ≔1⟨d(k),y(k)⟩⋅⟨y(k)−2‖y(k)‖2⟨d(k),y(k)⟩d(k),s(k+1)⟩ 1 d (k ),y (k ) · y (k )-2 \|y (k ) \|^2 d (k ),y (k ) d (k ),s (k+1 ) η(k) η (k ) ≔−1‖d(k)‖minη,‖s(k)‖ - 1 \|d (k ) \| \ \η, \|s (k ) \| \ β¯(k) β (k ) ≔maxβ(k),η(k) \ \β (k ),η (k ) \ It was proven that if the line search step in Algorithm 3 satisfies Wolfe-Powell conditions and the gradient is globally Lipschitz, Hager-Zhang conjugate gradient achieves global convergence finding a stationary point for a smooth nonconvex objective function. In a dynamical system view, this corresponds to the global attraction property of the trajectory of the numerical algorithm to find equilibria for globally Lipschitz flows. Lemma 11 implies that ρf~G_ρ f, or s(k)s (k ) in Algorithm 3, are globally Lipschitz. Thus, by Lemma 10, Algorithm 3 yields the Clarke stationary point of f. Based on the arguments in Section 4.2, if the potential function for ρf~G_ρ f is tractable to construct, the Wolfe-Powell line search in Algorithm 3 can be carried out using the potential function of ρf~G_ρ f as the surrogate objective function; alternatively, an exact line search can be carried out by finding α that satisfies ⟨ρf~(x(k)+α⋅d(k)),d(k)⟩=0 _ρ f (x (k )+α· d (k ) ),d (k ) =0 — such an exact line search can usually be carried out efficiently using Brent’s method to find a root of a one-dimensional equation in ℝ>0R_>0 7. Furthermore, the descent property of d(k)d (k ) was shown by 22 independent of the line searches, which guarantees that ⟨ρf~(x(k)+α⋅d(k)),d(k)⟩=0 _ρ f (x (k )+α· d (k ) ),d (k ) =0 has a positive root. Moreover, another line search to ensure global convergence can be carried out by backtracking to find α(k)α (k ) satisfying −⟨ρf~(x(k)+c1⋅α(k)d(k)),d(k)⟩≥c1c2⋅α(k)‖d(k)‖2,- _ρ f (x (k )+c_1·α (k )d (k ) ),d (k ) ≥ c_1c_2·α (k ) \|d (k ) \|^2, (92) where c1,c2∈ℝ>0c_1,c_2 _>0 are constant to be chosen. When ρf~G_ρ f is pseudo-monotone in the sense of Karamardian 30, since the global Lipschitz property was established for ρf~G_ρ f in Lemma 11, global convergence was proven for this backtracking line search method 16. We conclude this section with the observation that certain conjugate gradient methods obviate the need for line search procedures by determining the step size directly from s(k)s^(k) and d(k)d^(k), as exemplified in 8. 5 Optimizing Algorithm and Prediction for Penalized qqGaussian Likelihood Problems 5.1 Problem Formulation Using the qqGaussian distribution to model the data will undoubtedly enhance the robustness towards the underlying distributional assumption and outliers. However, unlike the Gaussian distribution, two independent qqGaussian random vectors are not jointly qqGaussian. Thus, we take the following approach to model the data. Let train∈ℝntrain×(p+1)X_train ^n_train× (p+1 ), ytrain∈ℝntrain\ y_train ^n_train denote the training design matrix and outcome, val∈ℝnval×(p+1),yval∈ℝnvalX_val ^n_val× (p+1 ),\ y_val ^n_val denote the validation design matrix and outcome, and test∈ℝntest×(p+1),ytest∈ℝntestX_test ^n_test× (p+1 ),\ y_test ^n_test denote the testing design matrix and outcome. Let ≔[trainT,valT,testT]T∈ℝn×(p+1)X [X_train^T,X_val^T,X_test^T ]^T ^n× (p+1 ) (93) denote the design matrix for the entire dataset, and let y≔[ytrainT,yvalT,ytestT]T∈ℝny [y_train^T,y_val^T,y_test^T ]^T ^n (94) denote the outcome for the entire dataset. Instead of assuming the qqGaussian distribution for the training, validation and testing set separately, we assume that y∼qGaussian(q,θ,Σ),y qGaussian (q,Xθ, ), (95) where θ∈ℝp+1θ ^p+1 denotes the coefficients for regression, and Σ denotes the characteristic/scale matrix for the entire data. Clearly, train _train =[Intrain×ntrain,0ntrain×nval,0ntrain×ntest] = [I_n_train× n_train,0_n_train× n_val,0_n_train× n_test ]X ytrain y_train =[Intrain×ntrain,0ntrain×nval,0ntrain×ntest]y = [I_n_train× n_train,0_n_train× n_val,0_n_train× n_test ]y implies that ytrain∼qGaussian(qtrain,trainθ,Σtrain)y_train qGaussian (q_train,X_trainθ, _train ) (96) where by the linear mapping closeness property 2, Σtrain=[Intrain×ntrain,0ntrain×nval,0ntrain×ntest]Σ[Intrain×ntrain,0ntrain×nval,0ntrain×ntest]T _train= [I_n_train× n_train,0_n_train× n_val,0_n_train× n_test ] [I_n_train× n_train,0_n_train× n_val,0_n_train× n_test ]^T (97) is the ntrain×ntrainn_train× n_train block diagonal matrix of Σ corresponding to the training data. By (2), 21−qtrain−1−ntrain=21−q−1−n, 21-q_train^-1-n_train= 21-q^-1-n, (98) which implies that 1qtrain−1−ntrain=1q−1−n. 1q_train-1-n_train= 1q-1-n. (99) (99) allows us to recover q from the training procedure. Above formulas for the training data and parameters can trivially be applied to the validation and the testing data and parameters; thus, validation and test can carried out easily from the model build from the training data. For q−q-correlated data, often times, the q−q-correlation structure is inferred or given prior to the model fitting; thus, we assume that the q−q-correlation structure is given as Ψ and we estimate the volatility / dispersion / scale parameter σ2>0σ^2>0 such that Σ=σ2Ψ. =σ^2 . (100) Trivially, Ψtrain _train is the block diagonal matrix of Ψ corresponding to the training data and Σtrain=σ2Ψtrain _train=σ^2 _train. We are now ready to formulate our likelihood loss function. To utilize qqGaussian distribution to model the q−q-correlated observations, we estimate the value of q such that q is allowed to vary, and the model will thus be more robust towards a wide class of distributions. Therefore, we choose to build the model using (20), since the dispersion matrix Λ (31) depends on q. We formulate our maximization of our log-likelihood function as the following from (96) and (20): dgroup* argmaxqtrain∈(1,1+2ntrain),θ∈ℝp+1,σ2∈ℝ>0log(1|σ2Ψtrain|1/2⋅Γ(1qtrain−1)Γ(1qtrain−1−ntrain2)⋅(2qtrain−1−ntrain)−ntrain2⋅(1+(2qtrain−1−ntrain)−1⋅⟨ytrain−trainθ,(σ2Ψtrain)−1(ytrain−trainθ)⟩)11−qtrain). q_train∈(1,1+ 2n_train),θ ^p+1,σ^2 _>0 \ ( 1 |σ^2 _train | 12· ( 1q_train-1 ) ( 1q_train-1- n_train2 )· ( 2q_train-1-n_train )^- n_train2· (1+ ( 2q_train-1-n_train )^-1· y_train-X_trainθ, (σ^2 _train )^-1 (y_train-X_trainθ ) ) 11-q_train). To address the high–dimensional data concerns, Oracle penalties are incorporated to carry out variable selection. To penalize the log-likelihood loss function to achieve variable selection, we formulate the following problem: dgroup* argminqtrain∈(1,1+2ntrain),θ∈ℝp+1,σ2∈ℝ>0−log(1|σ2Ψtrain|1/2⋅Γ(1qtrain−1)Γ(1qtrain−1−ntrain2)⋅(2qtrain−1−ntrain)−ntrain2⋅(1+(2qtrain−1−ntrain)−1⋅σ−2⋅(⟨ytrain−trainθ,Ψtrain−1(ytrain−trainθ)⟩+2ntrain∑j=2p+1w(θj)))11−qtrain) q_train∈(1,1+ 2n_train),θ ^p+1,σ^2 _>0 \ - ( 1 |σ^2 _train | 12· ( 1q_train-1 ) ( 1q_train-1- n_train2 )· ( 2q_train-1-n_train )^- n_train2 · (1+ ( 2q_train-1-n_train )^-1·σ^-2· ( y_train-X_trainθ, _train^-1 (y_train-X_trainθ ) +2n_train _j=2^p+1w ( _j ) ) ) 11-q_train) ⇔argminqtrain∈(1,1+2ntrain),θ∈ℝp+1,σ2∈ℝ>0n2logσ2−logΓ(1qtrain−1)+logΓ(1qtrain−1−ntrain2)+ntrain2log(2qtrain−1−ntrain)+1qtrain−1log(1+(2qtrain−1−ntrain)−1⋅σ−2⋅(⟨ytrain−trainθ,Ψtrain−1(ytrain−trainθ)⟩+2ntrain∑j=2p+1w(θj))) q_train∈(1,1+ 2n_train),θ ^p+1,σ^2 _>0 \ n2 σ^2- ( 1q_train-1 )+ ( 1q_train-1- n_train2 )+ n_train2 ( 2q_train-1-n_train ) + 1q_train-1 (1+ ( 2q_train-1-n_train )^-1·σ^-2· ( y_train-X_trainθ, _train^-1 (y_train-X_trainθ ) +2n_train _j=2^p+1w ( _j ) ) ) (101) In the above formulated problem, w is the Oracle penalty function, and we are not to penalize the intercept term. The 2ntrain2n_train multiplier is to ensure that the penalization effect is consistent with the number of training observations. We choose to put the penalty term together with the quadratic term without the variance scale parameter σ2σ^2 for two reasons: first, the optimization problem is more tractable under such problem formulation; second, we do not wish to let the value of σ2σ^2 perturb the degree of penalization. Comparing to penalized the log-likelihood directly, we choose to penalize the quadratic component directly as it is more tractable. It was shown that doing so will preserve Oracle properties 40 of penalized estimators. For the optimization procedure, we will proceed in a blockwise manner; i.e., we will optimize qtrain,θ,σ2q_train,θ,σ^2 separately in each iteration. More details will be given in the following subsections. 5.2 Minimizing with respect to qtrainq_train and σ2σ^2 With all the other parameters fixed, the sub-problem to minimize with respect to σ2σ^2 is argminσ2∈ℝ>0 σ^2 _>0 \ n2logσ2+1qtrain−1log(1+(2qtrain−1−ntrain)−1⋅σ−2CLOSE n2 σ^2+ 1q_train-1 (1+ ( 2q_train-1-n_train )^-1·σ^-2 (102) ⋅(⟨ytrain−trainθ,Ψtrain−1(ytrain−trainθ)⟩+2ntrain∑j=2p+1w(θj))) · ( y_train-X_trainθ, _train^-1 (y_train-X_trainθ ) +2n_train _j=2^p+1w ( _j ) )) which has a smooth objective function with respect to σ2σ^2. The first-order optimality condition n2= n2= 1qtrain−1 1q_train-1 ⋅(2qtrain−1−ntrain)−1⋅(⟨ytrain−trainθ,Ψtrain−1(ytrain−trainθ)⟩+2ntrain∑j=2p+1w(θj))σ2+(2qtrain−1−ntrain)−1⋅(⟨ytrain−trainθ,Ψtrain−1(ytrain−trainθ)⟩+2ntrain∑j=2p+1w(θj)) · ( 2q_train-1-n_train )^-1· ( y_train-X_trainθ, _train^-1 (y_train-X_trainθ ) +2n_train _j=2^p+1w ( _j ) )σ^2+ ( 2q_train-1-n_train )^-1· ( y_train-X_trainθ, _train^-1 (y_train-X_trainθ ) +2n_train _j=2^p+1w ( _j ) ) implies that the optimal value for the subproblem (102) takes minimizer σ2¯ σ^2 =(1qtrain−1/n2−1)⋅(2qtrain−1−ntrain)−1 = ( 1q_train-1/ n2-1 )· ( 2q_train-1-n_train )^-1 (103) ⋅(⟨ytrain−trainθ,Ψtrain−1(ytrain−trainθ)⟩+2ntrain∑j=2p+1w(θj))>0, · ( y_train-X_trainθ, _train^-1 (y_train-X_trainθ ) +2n_train _j=2^p+1w ( _j ) )>0, which is feasible. The feasible set for qtrainq_train is (1,1+2ntrain) (1,1+ 2n_train ), in this view, when ntrainn_train is large, the numerical stability will be an issue if minimization is carried out with respect to qtrainq_train directly. Thus, we choose to minimize with respect to 1qtrain−1∈(ntrain2,∞) 1q_train-1∈ ( n_train2,∞ ). First of all, we are to prove that such minimization is feasible. Lemma 12. The objective function (101) has a local minimizer in (ntrain2,∞) ( n_train2,∞ ) with respect to 1qtrain−1 1q_train-1. proof Since the objective function (101) is continuous and smooth with respect to 1qtrain−1 1q_train-1, we only need to analyze the derivative when 1qtrain−1↘0 1q_train-1 0 and 1qtrain−1→∞ 1q_train-1→∞. 1qtrain−1→∞ 1q_train-1→∞: Stirling’s formula states that limx→∞Γ(x)2πx(xe)x(1+O(x−1))=1. _x→∞ (x ) 2πx ( xe )^x (1+O (x^-1 ) )=1. (104) Thus, lim1qtrain−1→∞Γ(1qtrain−1)Γ(1qtrain−1−ntrain2)⋅(2qtrain−1−ntrain)−ntrain2=1 _ 1q_train-1→∞ ( 1q_train-1 ) ( 1q_train-1- n_train2 )· ( 2q_train-1-n_train )^- n_train2=1 (105) then lim1qtrain−1→∞−log(Γ(1qtrain−1)Γ(1qtrain−1−ntrain2)⋅(2qtrain−1−ntrain)−ntrain2)=0. _ 1q_train-1→∞- ( ( 1q_train-1 ) ( 1q_train-1- n_train2 )· ( 2q_train-1-n_train )^- n_train2 )=0. (106) We also have 1qtrain−1log(1+(2qtrain−1−ntrain)−1⋅σ−2CLOSE 1q_train-1 (1+ ( 2q_train-1-n_train )^-1·σ^-2 ⋅(⟨ytrain−trainθ,Ψtrain−1(ytrain−trainθ)⟩+2ntrain∑j=2p+1w(θj))) · ( y_train-X_trainθ, _train^-1 (y_train-X_trainθ ) +2n_train _j=2^p+1w ( _j ) )) =O((1qtrain−1)/log(1qtrain−1)), =O ( ( 1q_train-1 )/ ( 1q_train-1 ) ), which implies that this term will goes to infinity as 1qtrain−1→∞ 1q_train-1→∞. Thus, the objective function (101) goes to infinity as 1qtrain−1→∞ 1q_train-1→∞. 1qtrain−1↘ntrain2 1q_train-1 n_train2: Since Γ(1qtrain−1−ntrain2)→∞ ( 1q_train-1- n_train2 )→∞ as 1qtrain−1↘ntrain2 1q_train-1 n_train2. The penalized log-likelihood involving 1qtrain−1 1q_train-1 can be simplified as −logΓ(1qtrain−1)Γ(1qtrain−1−ntrain2)⋅((2qtrain−1−ntrain) - ( 1q_train-1 ) ( 1q_train-1- n_train2 )·( ( 2q_train-1-n_train ) OPEN+⟨ytrain−trainθ,(σ2Ψtrain)−1(ytrain−trainθ)⟩+2ntrain∑j=2p+1w(θj))11−qtrain→∞ + y_train-X_trainθ, (σ^2 _train )^-1 (y_train-X_trainθ ) +2n_train _j=2^p+1w ( _j )) 11-q_train→∞ as 1qtrain−1↘ntrain2 1q_train-1 n_train2. Thus, the subproblem to minimize with respect to 1qtrain−1 1q_train-1 is coercive on (ntrain2,∞) ( n_train2,∞ ). Coercivity implies that any minimizing sequence (1qtrain−1)j \ ( 1q_train-1 )_j \ must be contained within a bounded subset of (ntrain2,∞) ( n_train2,∞ ). Thus, Bolzano–Weierstrass theorem implies the existence of a convergent subsequence. Let (1qtrain−1)jk \ ( 1q_train-1 )_j_k \ be one such subsequence, and let (1qtrain−1)¯ ( 1q_train-1 ) be its limit. Since the subproblem has a continuous objective function with respect to 1qtrain−1 1q_train-1, the objective function is lower–semicontinuous and the value of the objective function at (1qtrain−1)¯ ( 1q_train-1 ) is less than or equal to the value of the objective function at (1qtrain−1)jk ( 1q_train-1 )_j_k for all k=1,2,…,∞k=1,2,…,∞. Thus, since (1qtrain−1)j \ ( 1q_train-1 )_j \ is a minimizing sequence, the value of the objective function at (1qtrain−1)¯ ( 1q_train-1 ) is less than or equal to the infimum of the objective function on (ntrain2,∞) ( n_train2,∞ ). Hence, since the entire minimizing sequence is contained in (ntrain2,∞) ( n_train2,∞ ), (1qtrain−1)¯∈(ntrain2,∞) ( 1q_train-1 )∈ ( n_train2,∞ ) solves the subproblem of minimizing with respect to 1qtrain−1 1q_train-1. Unlike σ2σ^2, the minimizer for 1qtrain−1 1q_train-1 is not in closed form, and the evaluation of the derivative with respect to 1qtrain−1 1q_train-1 can not be carried out efficiently. Thus, we apply Brent’s line-search method to optimize the 1qtrain−1 1q_train-1 subproblem 7. 5.3 Minimizing with respect to θ Minimizing with respect to θ involves a nonconvex smooth function and a convex nonsmooth function, which is termed a composite problem. In Section 4.3, we developed a proximal conjugate gradient algorithm for such composite optimization. In this part, we will establish an important remark regarding the Oracle penalty. The subproblem we are to minimize with respect to θ is argminθ∈ℝp+1⟨ytrain−trainθ,Ψtrain−1(ytrain−trainθ)⟩+2ntrain∑j=2p+1w(θj) θ ^p+1 \ y_train-X_trainθ, _train^-1 (y_train-X_trainθ ) +2n_train _j=2^p+1w ( _j ) ⇔ argminθ∈ℝp+112ntrain⟨ytrain−trainθ,Ψtrain−1(ytrain−trainθ)⟩+∑j=2p+1w(θj) θ ^p+1 \ 12n_train y_train-X_trainθ, _train^-1 (y_train-X_trainθ ) + _j=2^p+1w ( _j ) ⇔ argminθ∈ℝp+112ntrain⟨θ,trainTΨtrain−1trainθ⟩−2⟨ytrain,Ψtrain−1trainθ⟩+∑j=2p+1w(θj) θ ^p+1 \ 12n_train θ,X_train^T _train^-1X_trainθ -2 y_train, _train^-1X_trainθ + _j=2^p+1w ( _j ) (107) w can be chosen as oracle penalties such as SCAD/MCP penalties. And it has been shown that both SCAD / MCP penalties admit a difference-of-convex decomposition to a first-order smooth concave term plus λ times ℓ1 _1 penalty. The quadratic loss function is clearly convex and smooth. This justifies our assumption for the objective function. To carry out the proximal Hager-Zhang conjugate gradient method proposed in Section 4.3, we need to calculate L∇gL_∇ g, the L−L-smoothness constant for the smooth component. Previous work suggests L∇g=maxmax eigenvalue of 1ntraintrainTΨtrain−1train,cpenaltyL_∇ g= \ \max eigenvalue of 1n_trainX_train^T _train^-1X_train,c_penalty \, where cpenaltyc_penalty is the L−L-smoothness constant for the smooth component of the penalty, which will be 1a−1 1a-1 for SCAD and 1γ 1γ for MCP 63. Remark 13. For high dimensional data, often times, the number of covariates exceeds the number of observations; i.e, null(train)≠∅null (X_train )≠ . Both SCAD/MCP penalties take constant values in B∞(0,c)B_∞ (0,c )11 1 B∞(0,c)B_∞ (0,c ) denotes the open ball in uniform norm in the corresponding space, centered at the origin with radius c. ; where c=aλc=aλ for SCAD and c=γλc=γλ for MCP. Given any stationary point θ¯ θ in the nonempty solution set defined by trainTΨtrain−1train−trainTytrain=0X_train^T _train^-1X_train-X_train^Ty_train=0. For the set θ¯+null(train)∖B∞(0,c) θ+null (X_train ) B_∞ (0,c ), each point in the relative interior (which is nonempty) of this set is a Clarke stationary point, which implies that any algorithm with a starting point in this set will converge in 00 steps. This might pose an issue for signal recovery, since null(train)null (X_train ) is a vector subspace and some points can be very far from the origin. Remark 14. In view of the subproblem with respect to θ, it is trivial that the minimizer for (107) does not depend on the other parameters, which are q and σ2σ^2. Since the qqGaussian distribution is a generalization for all bell curve distributions, the estimation of the central trend using the maximum likelihood principle for bell curve distributions is equivalent to minimize a quadratic function, which has a breakdown point of 00. Taking into account the optimization subproblem with respect to θ, it is evident that the solution to (107) remains unaffected by the other parameters, namely q and σ2σ^2. Given that the qqGaussian distribution extends the framework of bell curve distributions, the problem (107) implies that estimating the central trend through the maximum likelihood principle for all bell curve distributions is equivalent to minimizing a quadratic function of the central trend, therefore characterized by a breakdown point of 00. 5.4 Prediction for ytesty_test To show how prediction can be made, we will show the methods to predict ytesty_test using the trained model in this subsection. The same method applies for validation when predictions on yvaly_val are needed or to predict any new data based on the trained model. Since the data are mutually qqGaussian, (99) implies that 1qtrain−1−ntrain=1qval−1−nval=1qtest−1−ntest=1q−1−n, 1q_train-1-n_train= 1q_val-1-n_val= 1q_test-1-n_test= 1q-1-n, (108) which will be used to recover the value of the shape parameter qval,qtestq_val,q_test. Note that when nnewn_new data points are introduced, the total number of observations n changes from ntrain+nval+ntestn_train+n_val+n_test to ntrain+nval+ntest+nnewn_train+n_val+n_test+n_new, thus, the value of q⋅q_· will change for the entire dataset. However, qtrainq_train stays the same; thus, we suggest inferring the shape parameter for each data set based on qtrainq_train directly using the equation above. With q⋅q_· calculated, it is straightforward to estimate the q−q-variance-covariance matrix q[(y⋅−θ)(y⋅−θ)T]E_q [ (y_·-Xθ ) (y_·-Xθ )^T ] based on (33); or, if existing, the variance-covariance matrix [(y⋅−θ)(y⋅−θ)T]E [ (y_·-Xθ ) (y_·-Xθ )^T ] based on (34). 6 Conclusion and Discussion This paper explores the field of statistical sparse learning, focusing on modeling correlated data through the lens of maximizing Tsallis entropy. It addresses the limitations inherent in the conventional Gaussian distribution, notably its lack of robustness towards outliers and underlying shape assumptions, by advocating for the qqGaussian distribution. This distribution, derived from Tsallis entropy maximization, represents a novel approach to handling correlated data and heterogeneity — elements frequently encountered in biostatistical contexts involving genetic and longitudinal studies. This paper encompasses a re-derived probability density function for the multivariate qqGaussian distribution based on Tsallis entropy maximization. Statistical modeling based on the derived density paves the way for the analysis of correlated data and heterogeneity and enables variable selection. Furthermore, we have developed an innovative framework capable of converting any numerical method, originally designed to identify equilibria in flows, into a tool for tackling composite optimization problems that are prevalent in statistical sparse learning. By applying this framework to the Hager-Zhang conjugate gradient algorithm, we have crafted an effective and stable algorithm tailored to the challenges of sparse statistical learning. Given the abundance of methods for numerically identifying equilibria for globally Lipschitz flows, our approach significantly broadens the arsenal of techniques available to address sparse statistical learning optimization challenges. In conclusion, our research positions the qqGaussian distribution, underpinned by maximizing Tsallis entropy, as a robust and adaptable alternative to Gaussian-based methodologies in statistical sparse learning on correlated data. This breakthrough not only confronts the traditional limitation of Gaussian assumptions, but also paves the way for expanded investigation into Tsallis entropy-maximizing distributions, particularly within the domain of biostatistics and allied disciplines. Future directions for research include the exploration of the log-linear model through the lens of Tsallis entropy maximization, akin to approaches previously based on Shannon’s entropy. Moreover, the study of the phenomenon called volatility smirk in financial return data may benefit from employing the log-qqGaussian distribution —- a transformation of the qqGaussian distribution, which can provide deeper insights into the nuances of financial markets. Additionally, in the field of statistical computing research, our framework that transforms numerical methods for identifying flow equilibria into algorithms for solving composite optimization problems opens numerous avenues for future research, especially in a sparse learning context. 7 Acknowledgments This work was supported by the ISM Scholarship for Outstanding PhD Candidates awarded to K. Yang, the NSERC Discovery Grant to C. Greenwood (Grant Number: RGPIN-2019-04482), the NSERC Discovery Grant to M. Asgharian (Grant Number: RGPIN-2024-05640), and the CANSSI Collaborative Research Team Grant to C. Greenwood and G. Cohen Freue. References Agarwal et al. (2009) R. P. Agarwal, M. Meehan, and D. O’Regan Fixed point theory and applications. Digitally printed version, paperpack re-issue edition, Cambridge tracts in mathematics, Cambridge Univ. Press, Cambridge [u.a.]. External Links: ISBN 9780521802505 Cited by: §4.2. Atkinson (1989) K. E. Atkinson An introduction to numerical analysis. 2. ed., [14. print] edition, Wiley, New York [u.a.]. Note: Bibliogr. S. 665 External Links: ISBN 0471624896 Cited by: §4.2, §4.2. Bauschke (2011) H. H. BauschkeP. L. Combettes (Ed.) Convex analysis and monotone operator theory in hilbert spaces. SpringerLink, Springer New York, New York, NY. External Links: ISBN 9781441994677 Cited by: §4. Beck and Teboulle (2009) A. Beck and M. Teboulle A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences 2 (1), p. 183–20 (English). Note: Copyright - Copyright] © 2009 Society for Industrial and Applied Mathematics; Last updated - 2012-07-02 External Links: Link Cited by: §4.2. Borland (2002a) L. Borland A theory of non-gaussian option pricing. Quantitative Finance 2 (6), p. 415–431. External Links: Document Cited by: §1. Borland (2002b) L. Borland Option pricing formulas based on a non-gaussian stock price model. Physical Review Letters 89 (9), p. 098701. External Links: Document Cited by: §1. Brent (1971) R. P. Brent An algorithm with guaranteed convergence for finding a zero of a function. The Computer Journal 14 (4), p. 422–425. External Links: ISSN 1460-2067, Document Cited by: §4.2, §4.3, §5.2. Chen et al. (2018) C. Chen, L. Luo, C. Han, and Y. Chen Global convergence of an extended descent algorithm without line search for unconstrained optimization. Journal of Applied Mathematics and Physics 06 (01), p. 130–137. External Links: ISSN 2327-4379, Document Cited by: §4.2, §4.3. Clarke (1990) F. H. Clarke Optimization and nonsmooth analysis. Classics in applied mathematics, Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104), Philadelphia, Pa. Note: Reprint. Originally published: New York : Wiley, 1983 External Links: ISBN 9781611971309 Cited by: §4. Clarke (2004) F. Clarke Lyapunov functions and feedback in nonlinear control. In Lecture Notes in Control and Information Sciences, p. 267–282. External Links: ISBN 9783540399834, Document, ISSN 1610-7411 Cited by: §4.2. Costa et al. (2003) J. Costa, A. Hero, and C. Vignat On solutions to multivariate maximum α-entropy problems. In Lecture Notes in Computer Science, p. 211–226. External Links: Document Cited by: §3, §3. Cover and Thomas (2006) T. M. Cover and J. A. Thomas Elements of information theory. John Wiley & Sons, Inc.. External Links: ISBN 9780471241959 Cited by: §1, §2. Dai and Yuan (1999) Y. H. Dai and Y. Yuan A nonlinear conjugate gradient method with a strong global convergence property. SIAM Journal on Optimization 10 (1), p. 177–182. External Links: ISSN 1095-7189, Document Cited by: §4.2, §4.3. Dandine-Roulland and Perdry (2015) C. Dandine-Roulland and H. Perdry The use of the linear mixed model in human genetics. Human Heredity 80 (4), p. 196–206. External Links: ISSN 1423-0062, Document Cited by: §1. Domingo et al. (2017) D. Domingo, A. d’Onofrio, and F. Flandoli Boundedness vs unboundedness of a noise linked to tsallis q-statistics: the role of the overdamped approximation. Journal of Mathematical Physics 58 (3). External Links: ISSN 1089-7658, Document Cited by: §1. Feng et al. (2017) D. Feng, M. Sun, and X. Wang A family of conjugate gradient methods for large-scale nonlinear equations. Journal of Inequalities and Applications 2017 (1). External Links: ISSN 1029-242X, Document Cited by: §4.2, §4.3. Fletcher (1964) R. Fletcher Function minimization by conjugate gradients. The Computer Journal 7 (2), p. 149–154. External Links: ISSN 1460-2067, Document Cited by: §4.2, §4.3. Garcia and Marder (2017) T. P. Garcia and K. Marder Statistical approaches to longitudinal data analysis in neurodegenerative diseases: huntington’s disease as a model. Current Neurology and Neuroscience Reports 17 (2). External Links: ISSN 1534-6293, Document Cited by: §1. Ghadimi and Lan (2013) S. Ghadimi and G. Lan Stochastic first- and zeroth-order methods for nonconvex stochastic programming. SIAM Journal on Optimization 23 (4), p. 2341–2368. External Links: Document Cited by: §4.2, §4.2. Ghadimi and Lan (2015) S. Ghadimi and G. Lan Accelerated gradient methods for nonconvex nonlinear and stochastic programming. Mathematical Programming 156 (1-2), p. 59–99. External Links: Document Cited by: §1. Gilbert and Nocedal (1992) J. C. Gilbert and J. Nocedal Global convergence properties of conjugate gradient methods for optimization. SIAM Journal on Optimization 2 (1), p. 21–42. External Links: ISSN 1095-7189, Document Cited by: §4.3. Hager and Zhang (2005) W. W. Hager and H. Zhang A new conjugate gradient method with guaranteed descent and an efficient line search. SIAM Journal on Optimization 16 (1), p. 170–192. External Links: ISSN 1095-7189, Document Cited by: item 3, §4.2, §4.3, §4.3, §4.3, §4.3, Abstract, Algorithm 3. Hager and Zhang (2006) W. Hager and H. Zhang A survey of nonlinear conjugate gradient method. 2. Cited by: §4.3. Helmke (1994) U. Helmke Optimization and dynamical systems. Springer London, London. External Links: ISBN 9781447134671, Document, ISSN 0178-5354 Cited by: §4.2. Hestenes and Stiefel (1952) M.R. Hestenes and E. Stiefel Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards 49 (6), p. 409. External Links: ISSN 0091-0635, Document Cited by: §4.2, §4.3. Hoheisel et al. (2020) T. Hoheisel, M. Laborde, and A. M. Oberman A regularization interpretation of the proximal point method for weakly convex functions. Journal of Dynamics & Games. External Links: Link Cited by: §1. Hubbard and West (1995) J. H. Hubbard and B. H. West Differential equations: a dynamical systems approach. Springer New York. External Links: ISBN 9781461241928, Document, ISSN 0939-2475 Cited by: §4.2. Kafka and Wilke (2019) D. Kafka and D. Wilke Gradient-only line searches: an alternative to probabilistic line searches. External Links: Document, 1903.09383 Cited by: §4.2. Kanzow and Lechner (2020) C. Kanzow and T. Lechner Globalized inexact proximal newton-type methods for nonconvex composite functions. Computational Optimization and Applications 78 (2), p. 377–410. External Links: ISSN 1573-2894, Document Cited by: §4.2. Karamardian (1976) S. Karamardian Complementarity problems over cones with monotone and pseudomonotone maps. Journal of Optimization Theory and Applications 18 (4), p. 445–454. External Links: ISSN 1573-2878, Document Cited by: §4.3. Lee et al. (2014) J. D. Lee, Y. Sun, and M. A. Saunders Proximal newton-type methods for minimizing composite functions. SIAM Journal on Optimization 24 (3), p. 1420–1443. External Links: Document Cited by: Remark 5. Li et al. (2016) X. Li, H. Jiang, J. Haupt, R. Arora, H. Liu, M. Hong, and T. Zhao On fast convergence of proximal algorithms for sqrt-lasso optimization: don’t worry about its nonsmooth loss function. arXiv. External Links: Document Cited by: §4.2. Lubich (2006) C. LubichG. Wanner and E. Hairer (Eds.) Geometric numerical integration. 2nd ed. edition, Springer Series in Computational Mathematics Ser., Springer Berlin / Heidelberg, Berlin, Heidelberg. Note: Description based on publisher supplied metadata and other sources. External Links: ISBN 9783540306665 Cited by: §4.2. M. Tsukada (2005) M. K. M. Tsukada On the probability distribution maximizing generalized entropies. In Proceedings of 2005 Symposium on Applied Functional Analysis - Information Sciences and Related Fields, p. 99–111. Cited by: §3, §3. Morduchovič (2018) B. S. Morduchovič Variational analysis and applications. Softcover re-print of the Hardcover 1st edition 2018 edition, Springer monographs in mathematics, Springer, Cham, Switzerland. Note: Literaturverzeichnis: Seiten 533-578 External Links: ISBN 9783030065133 Cited by: §4. Mordukhovich (2006a) B. S. Mordukhovich Variational analysis and generalized differentiation i: basic theory. Variational analysis and generalized differentiation, Springer, Berlin ;. Note: Includes bibliographical references and indexes. External Links: ISBN 9783540312475 Cited by: §4. Mordukhovich (2006b) B. S. Mordukhovich Variational analysis and generalized differentiation i: applications. Variational analysis and differentiation, Springer, New York. Note: Includes bibliographical references and index. External Links: ISBN 9783540312468 Cited by: §4. Nesterov (2004) Y. Nesterov Introductory lectures on convex optimization. Springer US. External Links: Document Cited by: §4.2. Neumaier et al. (2024) A. Neumaier, M. Kimiaei, and B. Azmi Globally linearly convergent nonlinear conjugate gradients without wolfe line search. Numerical Algorithms. External Links: ISSN 1572-9265, Document Cited by: §4.2. Nikolova (2000) M. Nikolova Local strong homogeneity of a regularized estimator. SIAM Journal on Applied Mathematics 61 (2), p. 633–658. External Links: Document Cited by: §1, §4.2, §5.1. Nocedal and Wright (2000) J. Nocedal and S. WrightS. J. Wright (Ed.) Numerical optimization. Second edition edition, Springer Series in Operations Research and Financial Engineering, Springer New York, New York, NY. External Links: ISBN 9780387987934, LCCN 99013263, Link Cited by: §1. Peña et al. (1999) I. Peña, G. Rubio, and G. Serna Why do we smile? on the determinants of the implied volatility function. Journal of Banking & Finance 23 (8), p. 1151–1179. External Links: ISSN 0378-4266, Document Cited by: §1. Polak and Ribiere (1969) E. Polak and G. Ribiere Note sur la convergence de méthodes de directions conjuguées. ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique 3 (R1), p. 35–43 (fre). External Links: Link, MathReview Entry Cited by: §4.2, §4.3. Quarteroni et al. (2007) A. Quarteroni, R. Sacco, and F. Saleri Numerical mathematics. Springer New York. External Links: ISBN 9780387227504, Document, ISSN 2196-9949 Cited by: §4.2. Riahi and Qattan (2018) M. K. Riahi and I. A. Qattan Linearly convergent nonlinear conjugate gradient methods for a parameter identification problems. External Links: Document, 1806.10197 Cited by: §4.2, §4.2. Rockafellar and Wets (2010) R. T. Rockafellar and R. J.-B. Wets Variational analysis. Corr. 3. printing. [Softcover version of original hardcover edition 1998] edition, Die @Grundlehren der mathematischen Wissenschaften in Einzeldarstellungen, Springer, Heidelberg. External Links: ISBN 3642083048 Cited by: §4. Ross (2019) I. M. Ross An optimal control theory for accelerated optimization. External Links: Document, 1902.09004 Cited by: §4.2, §4.2. Runcie and Crawford (2019) D. E. Runcie and L. Crawford Fast and flexible linear mixed models for genome-wide genetics. PLOS Genetics 15 (2), p. e1007978. External Links: ISSN 1553-7404, Document Cited by: §1. Saad (2003) Y. Saad Iterative methods for sparse linear systems. Society for Industrial and Applied Mathematics. External Links: ISBN 9780898718003, Document Cited by: §4.2. Shannon (1948) C. E. Shannon A mathematical theory of communication. Bell System Technical Journal 27 (3), p. 379–423. External Links: Document Cited by: §2, §2. Shi and Shen (2005) Z. Shi and J. Shen Convergence of descent method without line search. Applied Mathematics and Computation 167 (1), p. 94–107. External Links: ISSN 0096-3003, Document Cited by: §4.2. Snyman (1985) J. A. Snyman UNCONSTRAINED minimization by combining the dynamic and conjugate gradient methods. Quaestiones Mathematicae 8 (1), p. 33–42. External Links: ISSN 1727-933X, Document Cited by: §4.2. Snyman (2004) J. A. Snyman A gradient-only line search method for the conjugate gradient method applied to constrained optimization problems with severe noise in the objective function. International Journal for Numerical Methods in Engineering 62 (1), p. 72–82. External Links: ISSN 1097-0207, Document Cited by: §4.2. Sontag (1998) E. D. Sontag Mathematical control theory. Second Edition edition, Springer eBook Collection, Springer New York, New York, NY. External Links: ISBN 9781461205777, Document, ISSN 2196-9949 Cited by: §4.2. Sun and Zhang (2001) J. Sun and J. Zhang Global convergence of conjugate gradient methods without line search. Annals of Operations Research 103 (1/4), p. 161–173. External Links: ISSN 0254-5330, Document Cited by: §4.2. Tsallis (1988) C. Tsallis Possible generalization of boltzmann-gibbs statistics. Journal of Statistical Physics 52 (1-2), p. 479–487. External Links: Document Cited by: §2, §2. Vignat et al. (2004) C. Vignat, A. Hero I, and J. Costa About closedness by convolution of the tsallis maximizers. Physica A: Statistical Mechanics and its Applications 340 (1-3), p. 147–152. External Links: Document Cited by: §3, §3, §3. Vignat and Plastino (2005) C. Vignat and A. Plastino The p-sphere and the geometric substratum of power-law probability distributions. Physics Letters A 343 (6), p. 411–416. External Links: Document Cited by: §3. Vignat and Plastino (2007) C. Vignat and A. Plastino Scale invariance and related properties of q-gaussian systems. Physics Letters A 365 (5-6), p. 370–375. External Links: Document Cited by: §3. Vignat and Plastino (2009) C. Vignat and A. Plastino Why is the detection of q−q-gaussian behavior such a common occurrence?. Physica A: Statistical Mechanics and its Applications 388 (5), p. 601–608. External Links: ISSN 0378-4371, Document Cited by: §3. Wang (2006) C. Wang Some remarks on conjugate gradient methods without line search. Applied Mathematics and Computation 181 (1), p. 370–379. External Links: ISSN 0096-3003, Document Cited by: §4.2. Wu (2011) Q. Wu A nonlinear conjugate gradient method without line search and its global convergence. In 2011 International Conference on Computational and Information Sciences, External Links: Document Cited by: §4.2. Yang et al. (2024) K. Yang, M. Asgharian, and S. Bhatnagar Accelerated gradient methods for sparse statistical learning with nonconvex penalties. Statistics and Computing 34 (1). External Links: ISSN 1573-1375, Document Cited by: §1, §4.2, §5.3. Yu and Peng (2017) Y. Yu and J. Peng The moreau envelope based efficient first-order methods for sparse recovery. Journal of Computational and Applied Mathematics 322, p. 109–128. External Links: ISSN 0377-0427, Document Cited by: §4.2. Zhou (2009) G. Zhou A descent algorithm without line search for unconstrained optimization. Applied Mathematics and Computation 215 (7), p. 2528–2533. External Links: ISSN 0096-3003, Document Cited by: §4.2.