Weighted L1 and L0 Regularization Using Proximal Operator Splitting Methods
Authors: Zewude A. Berkessa
Patrik Waldmann
Venue: Accepted by TMLR ()

Abstract: This paper develops a joint weighted $\displaystyle \normlone$- and $\displaystyle \normlzero$-norm (WL1L0) regularization method by leveraging proximal operators and translation mapping techniques to mitigate the bias introduced by the $\displaystyle \normlone$-norm in applications to high-dimensional data. A weighting parameter $\alpha$ is incorporated to control the influence of both regularizers. Our broadly applicable model is nonconvex and nonsmooth, but we show convergence for the alternating direction method of multipliers (ADMM) and the strictly contractive Peaceman–Rachford splitting method (SCPRSM). Moreover, we evaluate the effectiveness of our model on both simulated and real high-dimensional genomic datasets by comparing with adaptive versions of the least absolute shrinkage and selection operator (LASSO), elastic net (EN), smoothly clipped absolute deviation (SCAD) and minimax concave penalty (MCP). The results show that WL1L0 outperforms the LASSO, EN, SCAD and MCP by consistently achieving the lowest mean squared error (MSE) across all datasets, indicating its superior ability to handling large high-dimensional data. Furthermore, the WL1L0-SCPRSM also achieves the sparsest solution.

Published in Transactions on Machine Learning Research (11/2024)
Weighted L and L Regularization Using Proximal Operator Splitting Methods
Zewude A. Berkessa
zewude, wwwudee.ereesssuoouu..
Research Unit of Mathematical Sciences University of Oulu
Patrik Waldmann
potrik. marrik.ualdmonnBBouu.
Reseorch Unit of Mathematical Sciences University of Oulu
Reviewed on OpenReview: https: hptp///pppee/eeeee//////e/ net/ forum? id= 9m2k96cDMK
Abstract
This paper develops joint weighted L'- and -nrm (WL1L0) regularization method by leveraging proximal operators and translation mapping techniques to mitigate the bias introduced by the Lnorm in applications to high-dimensional data. A weighting parameter a is incorporated to control the influence of both regularizers. Our broadly applicable model nonconvex and nonsmooth, but we show convergence for the alternating direction method of multipliers (ADMM) and the strictly contractive (SCPRSM). Moreover, we evaluate the effectiveness of our model both simulated and real high-dimensional genomic datasets by comparing with adaptive versions of the least absolute shrinkage and selection operator (LASSO), elastic net (EN), smoothly clipped absolute deviation (SCAD) and minimax concave penalty (MCP). The results show that WL1LO outperforms the LASSO, EN, SCAD and MCP by consistently achieving the lowest mean squared error (MSE) across all datasets, indicating its superior ability to handling large high-dimensional data. Furthermore, the WLILO-SCPRM also achieves the sparsest solution. Julia code for the WLILO-ADMM and WLILO-SCPRSM available https: //g and--(PRS))
Airceman-Rachfort splitting method
1 Introduction
High-dimensional statisties is rapidly growing field of research that focuses on statistical analysis in the presence of large number of variables predictors (p), often much larger than the sample size (n). For example, high-throughput measurements in genomics contain thousands or millions of variables, such single mucleotide polymorphism (SNP) markers and gene expression data for each individual. In such settings, traditional statistical methods often fail due to issues like overfitting, multicollnnearity and computational complexity. In recent years, a number of regularization mothods have been developed that impose a penalty on the size of the regression coefficients, which encourages sparsity and reduces the number of variables the model (Fan al., 2011 Fan & Lv 2010: Heinze al. 2018 Sparse learning techniques are essential analyzing high-dimensional data for increased prediction accuracy, reduced computational complexity and
enhanced interpretability of the results (Bühlmann & Van De Geer 2011 Giraud 2015
Wainwright 2019
Among various sparsity-inducing methods, the L' regularizer (also known as the least absolute shrinkage and selection operator (LASSO)) stands out for its convex nature and computational efficiency [Tibshirani 1996). It adds a penalty term to the loss function proportional to the absolute value of the regression coefficients normm,,, which tends to shrink the coefficients towards zero and force some coefficients exactly zero. Shrinking these coefficients helps avoid overfitting, which can happen when model memorizes the training data too well and does not perform well data. The LASSSO improves the accuracy of predictions

---

Published in Transactions on Machine Learning Research (11/2024)
unseen data by simultancously selecting features and mitigating overfitting. However, when coefficients are being shrunk, bias is introduced due to the bias-variance trade-off that is unavoidable in statistical learning. Furthermore, the LASSO tends to favor keeping larger coefficients over smaller ones which can lead to bias towards larger cocfficients in the model estimation process (Hastic et al. 2015).
In the specific context of genomic data, where the goal is often to identify genes associated with certain traits or diseases, this inaccurate selection can lead to the inclusion of incorrect genes in the model. This also poses risk in terms of impaired prediction, as the estimated coefficients of the selected genes contribute to predicting the trait of interest (Fan et al., 2014b Fan & Li, 2001 Johnstone & Titterington 2009 Toloşi & Lengauer 2011). Hence, the LASSO requires the fulfilment of the irrepresentable condition to obtain valid estimations (Zhao & Yu 2006). In cases where the underlying datasets fail to meet this condition, the LASSO method may not accurately select the appropriate variables, leading to incorrect discoveries and wrong conclusions. In practice, implementing the irrepresentable condition can be challenging. Studies show that nonconvex regularizers such as SCAD and MCP reduce bias and have better prediction propertics than Zhang 2010
the L,- regularizer (Bertsimas et al. 2020 Fan& & Li
2001
On the other hand, L" regularization, which is also known best subset selection (Hocking & Leslie, 1967) directly penalizes the number of non-zero coefficients in the model. It encourages sparsity, meaning it tends to produce models with fewer non-zero coefficients without any shrinkage. This results in a model that only includes the most relevant variables, simplifying the model and potentially improving its predictive performance by reducing overfitting. However, finding the optimal subset of variables using the L -o-m is computationally expensive because the L' is nonconvex. While L' regularization is commonly used because of its convex nature, the -n"norm is computationally expensive and often intractable, and hence not frequently used on large data sets (Hastie et al. 2020).
Another regularization method is L2, also known as ridge regularization, which shrinks the coefficients towards zero without eliminating any of them completely (Hoerl & Kennard 1970). Unlike the LASSO, ridge regression produces dense estimated regression coefficients, which means it does not perform feature selection. It reduces the size of the coefficients but does not. drive any of them to exactly zero. Furthermore, the elastic net (EN) is another regularization technique that combines the properties of ridge regression and LASSO regression (Zou& & Hastic 2005). It is particularly useful for datasets with many features, especially when some are highly correlated. The LASSO may struggle with grouped variable selection, often picking just one from correlated variables, whereas EN improves feature selection in such cases (Hastie et al. 2015
In this paper, we propose combining L and L regularization into a method denoted WL1LO for improved prediction and variable selection. I," regularivation encourages sparsity by shrinking some coefficients to zero, which helps reduce overfitting and is computationally efficient due to its convex properties. On the other hand, 7 regularization enforces strict sparsity, directly penalizing the number of non-zero coefficients (i.e., eliminating variables with negligible impact), leading to more interpretable models. This synergy offers a better balance between interpretability and predictive accuracy, particularly in high-dimensional settings like genomics. Furthermore, since L is an unbiased estimator and L' often introduces biases in estimation, L can be regarded as debiasing L in this setting, We achieve this goal by using a common regularization parameter and introducing a weight parameter that balances the importance of the two regularization methods.
We address the computational challenges that arise from optimizing the L norm by using proximal split- ting methods, translation mapping and the efficient optimization algorithms ADMM and SCPRSM. The prediction and sparsity properties of the WL1LO method is evaluated on one simulated and two real genomic datasots and compared with the popular LASSO, EN, SCAD and MCP regularizers.
2 Related Work
In the rapidly evolving landscape of technology and data, prediction has become a cornerstone for making informed decisions across various domains. Regularization techniques are pivotal in enhancing the per- formance and generalizability of predictive models, particularly when dealing with complex datasets and high-dimensional data. By imposing penalties on the model parameters, regularization helps prevent over-
2

---

Published in Transactions on Machine Learning Research (11/2024)
fitting, ensuring that the model captures the underlying patterns in the data. In this section, we will review key related works on regularization methods, highlighting significant advancements and methodologies.
We start by introducing standard regression model
y Xb+c.
(1)
where y IR" is the response vector, X E Rxx is the predictor matrix, b E IP is the vector of regression coefficients, and e E IR" is a noise (error) vector. For a vector b, we write the q-norm notation as
if = 0, Σ, 1 (b, 0),
||6||
C 16/10//00
ir < < <0,
max; [byl. if = oc,
where i 1, ... Here, the | is the -orm that is the number of nonzero elements in b. It is noteworthy that the -nom does not meet the criteria of a norm, specifically lacking the homogeneity property (Bock) 2017). Despite this, the term is widely used the literature, and for the sake of consistency, will retain its adoption.
Ridge regression, also known Tikhonov regularization was introduced by Hoerl & Kennard 1970 uses an L penalty term that shrinks all the coefficients and reduces their magnitudes. The ridge regression can be formulated as
argmin y Xb112 + +1
(2)
where À > 0 is regulazation parameter need to be tuned. This method is particularly effective in addressing multicollinearit in linear regression models. However, it does not necessarily set any coefficients to zero. Hence, ridge regression does not produce a sparse solution of estimated coefficients. On the other hand, LASSO regression
b argmin y xb112 A|
(3)
incorporates an L' regularization penalty, which encourages sparsity in the solution by setting some coeffi- cients exactly to zero (Tibshirani 1996).
Other penalty functions are introduced to provide a balance between inducing sparsity and reducing estima- tion bias, aiming to solve the optimization problem as
argmin by xb + Px
For example, the smoothly clipped absolute deviation (SCAD) penalty function was introduced by Fan & Li (2001) as an improvement over LASSO regularization, particularly for bias reduction. The SCAD penalty function defined as
1a b
if |b| S X, bl < ad
PSCAD
if
(5)
2(a-11
11+1)A
6 aA,
where À > 0 and 0 are unknown parameters. for various problems, and À needs to be tuned.
Fan & Li
2001 suggested that 3.7 is good choice
The minimax concave penalty (MCP) is another lype of penalty function introduced by Zhang 2010 The MCP penalty function defined as
62
if < Åa,
PMP (b)
2e
(6)
axa
b > aX.
2
According the estimation theorems of Zhang (2010), 3 is a good choice for MCP, and À still needs be tuned. MCP was developed to address the estimation bias of the LASSO and is generally easier optimize computationally compared to SCAD.
3

---

Published in Transactions on Machine Learning Research (11/2024)
Both SCAD and MCP aim to eliminate unimportant variables while preserving important ones, achieving the 'oracle property' as the sample size grows (n + oo). They both asymptotically select the correct model and produce normal, accurate coefficient estimates. MCP is effective with many sparse predictor groups but struggles with tightly clustered non-zero coefficients while SCAD has weaker grouping behavior compared to MCP (Ogutu Piepho 2014p. We maintain the use of 3.7 for SCAD and for MCP throughout the paper.
For the L regularization (best subset selection (BSS)), (b) can be written as
pBss (b) = Σ 1 (b, #0).
Exact optimization of problem with the I,0nnorm, as defined in 7 is challenging because incorporating into the objective function results in a non-differentiable and non-conver problem. For example, Louizos et al. (2017) propose a method for optimizing a relaxed version of the L norm for parametric models using a distribution called the hard concrete distribution (Maddison et al. 2016). which facilitates gradient- based optimization.
Yun et al. (019)) use family of M-estimators with trimmed regularization for general high-dimensional problems. The trimmed regularization problem can be formulated for LASSO as
6, A argmin y b,
Xb11 +
R lb:
subject
1T2p-1 mm2p-1 p h
[0, 1].
where h denotes the trimming parameter, which must be appropriately tuned, for instance, through cross- validation.
Another example is the elastic net (EN) regression which combines both L and L regularization penalties, providing a balanced approach to prediction accuracy on future data and model interpretation in linear regression models. It is formulated as
6 argmin byy Xb12 + A1|16|I + 12116122
(9)
which has two regularization parameters X1 and A2 to tune [oo & Hastie 2005 The LAVA regression
model based on the splitting of the regression component into one sparse and one dense part b = c - d and thereby obtaining the following optimization problem
é, = arrmmnn y - X(ce d) + Alllell1 + +2 c.d
(10)
where the resulting estimator b d (Chenoozhukov et al. 2017 ). The key difference between EN and LAVA is that EN performs variable selection (i.c., is dominated by the L'-norm), whereas LAVA is always dense (i.e., is dominated by the n.omn. Waldmann (2021) developed a proximal operator algorithm based on the LAVA regularization method that jointly performs L.'- and --2-norm regularization.
Ziyin & Wang 2023 propose a method called spred, for optimizing generic differentiable objectives with an L constraint using a reparametrization. The method is proposed to effectively bridge the gap between sparsity in deep learning and conventional statistical learning by providing a principled way to optimize L constraints in complex nonlinear settings. For example, one can apply the spred parametrization to the | Xbb113 + 2x|[bbll1. The equivalent spred loss is then
sparse component of b, given the LASSO loss
é, d = argmin by ( 1) + K< [|] + llcl
(11)
e.d
where b= co å, O denotes the element-wise product, and ε is the L' regularization strength parameter that needs to be tuned achieve the best maassty-perforranne trade-off.
4

---

Published in Transactions on Machine Learning Research (11/2024)
3 Theoretical Background
Large parts of the theory behind our approach follows from (Bertsekas 2016 and
BBec
2017
For
extended real-valued function f RP + |-o0, oo], we define the following:
(a) The domain of f the set
dom(f) = (b R f(b) < co).
(b) is proper if dom(f) Ø and is never -o0.
(c) The epigraph of f is defined by
epi(f) ((b,a) RD x R: f(b) a).
(d) The function f is closed if its epigraph is closed.
(e) f is called lower semicontinuous at E R if
((b) S lim ((b(k)
for any sequence (b()) 1221 C R2 for which ((k) - b as h
(f) For any 7) E R, the n-level set of function f is the sct
Lev(J,n) (b R f(b) $n)
(g) A proper function f is called coercive if
lim ||b|]~000
f(b)
For any set 8 C R2 and any point b E RP, the distance from b to S is defined as D(b,8) = inf 1 bl1. m e 8}, and (bb, S) = oo for all when S = Ø.
A proper closed and cocrcive function f attains its minimal value over for a nonempty closed set satisfying sin dom(/) 0. Moreover, closed coercive function possesses a minimizer on any closed set that has a
nonempty intersection with the domain of the function (Beck) 2017).. For extended real-valued function f R -o s0|, the following three claims equivalent:
f is lower semicontinuous.
ii closed.
iii For any n E F, the level set
Lev(f.n) = (b E RP f(b) s n)
closed.
The proof of these claims can be found in (Beck 2017
Theorem 2.6.
3.1 Subdifferentials of Nonconvex and Nonsmooth Functions
Subdifferentials are important in analyzing complex functions, especially when dealing with nonsmooth and nonconvex functions. Following Clarke al. (2008) Mordukkoooi (2006), we explore subdifferentiabiiity.
Let g: RPP (-o0, to0] be proper and lower semicontinuous function. Then

---

Published in Transactions on Machine Learning Research (11/2024)
(i) For given b dom 9, the Fréchet subdifferential of at b, denoted by åg(b). is the set of all vectors IR? which satisly
lim inf g(m) g(b) (u, m b) 2
m-
m bl
and we set dg = when b $ dom g.
(ii The limiting-suubr saantin-sabdiffeeenntial, or simply the subdifferential, of g at b, written by g(b), is defined by
0g(b) = {u E R (b(k) - b, b(b((()) - g(b) and (h E @g(b(k)) k= k3 u),
where (b( c (b) for all IRP.
(ii)) A point b* is called critical point stationary point of if it satisfies 0 E 0g(b").
Please refer to Wu et al. [2021] for generalized subdiffercntials the L with its rogular subdiffercntials provided
Le
(2013).
3.2 The Kurdyka-tojasicn
Inequality and its Property
(KL) inequality deals with the behavior of certain functions near their critical points. It is an important tool for analyzing the convergence of nonconvex nonsmooth optimization prob- lems (Attouch et al. 2010 2013 Bolte al. 2014). We now revice the Kt property.
The
Let g: R2 (-o0, too be a proper lower semicontinuous function. Then,
(a) The function g Rº + R U (+o0} is said to have the KE property at b* E dom Oy if there exist E +oo], a neighborhood U of b", and contimuous concave function o [0. m)) - R such that
(i) (00) 0,
(ii) & is continuously disferentiable on (0,
(iii) Va (0, +oc], ('(a)
(IV) For all b U (g(b') < g(b) < g(b") + m} the KE property holds: (g((b) g )d((0 0g(b)) 2 1.
(b) Proper lower semicontinuous functions which satisfy the KL inequality at each point of dom åg are called KL functions. Examples of KE functions include |[b|],, |blloo and ly X For more examples, please refer to (Attouch et al.
2010 2013 Bolte al. 2014 Yashtini 2022
3.3 Proximal Operators
Proximal operators a fundamental concept in optimization, especially for problems involving non- smooth or non-convex functions, which are increasingly common in wide range of real-world applica- tions (Fukushima & Mine 1981 Kaplan & Tichatschke 1998 Parikh & Boyd 2013). proximal operator, denoted as prox, (u), aims to find a point closer to u that also minimizes specific objective function, f(v)
in specific optimization subproblem. This subproblem assumed to be more manageable to solve than the original problem. The proximal operator can be mathematically expressed as
proxy(1) = argmin (f(v) + (1/2)]l - -
(12)
where and are vectors of length Here, prox, (u) is point that compromises between minimizing
f and being close to u. Note that the right-hand side of (12) is strongly convex, hence there unique minimizer for every u E RR.. Introducing the parameter 7 > 0 that represents a trade-off parameter between
1
1
the two terms v and u yields scaled version of (12), in which
is replaced by 2-
The proximal operator
2
has useful properties [Bck 2017). one of which behavior when applied to affine functions, as shown below.

---

Published in Transactions on Machine Learning Research (11/2024)
Lemma 1 For any offine function f(u) (m, u) + a, where RR is fixed vector and E R, then for any E IR?, the prozimal operator defined in 12) reduces to simple translation of the veclor u by m. Specifically,
( = u - m,
(13)
which represents a translation mapping.
The proof of Lemma is provided in Appendix B. In accordance with Lemma one defines a translation
function as a function that incorporates a standard additive term which is expressed as (m) ) = m)
m. Another important property arises in the context of separable sum functions [(u, = ((u) + h(m),
where the proximal operator written as prox, m) = prox, + prox, (m). For proximal operators in
the framework of L", please refer to (Attouch et al. 2013 Bolte et al. 2014
Beck 2017).
4 Methodological Framework
To achieve higher sparsity than the EN, one can use the least squares loss function with é and to norm constraints which can be formulated as
b = argmin by Xb s.t. ||b|] S t, ||b|lo S
(14)
where 6 represents the estimate of the vector of regression coefficients, I is a constant threshold and is the desired level of sparsity (i.e., the maximum number nonzero coefficients). The optimization problem (14 can be expressed in the Lagrangian form as
= arggmn llainllllxett - Xb12 + 1116b11 + A2|16|lo.
(15)
where A > and X2 > 0 are the regularization parameters.
We formulate our proposed method by building upon the flexible, penalty-based framework introduced 15) that uses the parameters X and A2 to control both the size of the coefficients and the sparsity in a more flexible and nuanced manner, allowing for broader range of model behaviors compared the constant threshold and the strict subset selection enforced by formulation [4h. Here, the exact relationship between and λ₁ and between and X2 is data-dependent.
We now extend 15 by introducing a weight parameter E (0, 1), common regularization parameter À > 0 and reformulate the problem as
b = argmin by Xb + (aallbl + (1-a) (blob)
(16)
b
Furthermore, formulate the problem 16 following the LAVA method. Hence, we split the regression component b into the sparse components and d and separately assign penalties L and L to them and therefore obtain
WL1LO: è, argminn c.d
argmin 1y 1 X(c d)(l + X (allc]]) + (1-a)||d|))
(17)
It is here useful to point out that EN combines T.'-norm and T.2-ormm regularization which leads to variable selection while being less sensitive correlated predictors than the LASSO. In contrast, LAVA is dominated by the 1,--norm, leading to dense models without feature selection. WL1LO combines L and L regu- larization, providing stricter feature selection by explicitly controlling the number of non-zero coefficients, resulting in more precise sparsity than the EN. See Appendix A for additional insights into the problem and discussion.
5 Optimization Algorithms
Alternating direction method of multipliers (ADMM) is often used as a benchmark algorithm for splitting
problems due to its efficiency and flexibility (Boyd al. 2011 One of its notable strengths is its ability to

---

Published in Transactions on Machine Learning Research (11/2024)
handle large-scale optimization problems by decomposing them into smaller, more manageable subproblems. This decompasition not only simplifies the problem-solving process but also allows for parallel processing of these subproblems. The strictly contractive Peaceman Rachford splitting method (SCPRSM) is a variant of the classical Geaceman-Rachfone splitting method (PRSM) (He et al.) 2014). Similar to ADMM, PRSM is an operator splitting technique used solve optimization problems. SCPRSM a further exlension that ensures convergence by imposing a strict contraction condition. Although ADMM is more widely used practico, SCPRSM is more specialized method that guarantees convergence through striet contraction. While SCPRSM emphasizes strict contraction, ADMM may exhibit slower convergence under certain conditions. Therefore, we implement our proposed method using both ADMM and SCPRSM frameworks based on the augmented Lagrangian method that combines the original objective function with the constraints of the optimization problem into a single function. Here, the augmented Lagrangian's advantage lies in enabling the study of convergence for the proposed methods without requiring assumptions like strict convexity (Boyd ct al. 2011
The introduction of two variables (c and d) instead of one (b) in 17 increases the dimensionality of the optimization problem, adding complexity to the theoretical analysis. Therefore, we first study the convergence properties of the problem in (16). We then reformulate problem (16) to establish its convergence
in the ADMM and SCPRSM frameworks. The optimization model for [6 can be formulated as
b = argmin {f(b) + ((b)}
argmin (f(b) + g())
(18)
subject to b b= u,
where f(b) 11 x(b) is the loss function and g(b) (a |b|1 (1 -a)lbbllo)) the penalty function.
We now write the augmented Lagrangian function corresponding to 18
(b u, 2) s(b) + ((u) + =T u) + |1b u 1|
(19)
where z is a dual variable or Lagrange multiplier and > 0 is learning rate. Here, b and u are the primal variables.
5.1 Method of Multipliers and ADMM Framework
The method of multipliers jointly minimizes the two primal variables whereas the ADMM efficiently solves optimization problems by alternately updating primal and dual variables, elfectively decomposing complex problems into manageable subproblems (Boyd et al.) 2011). A more convenient scaled form of (19) can be obtained by completing the square with the dual variable and the residual b in the augmented Lagrangian. This allows the term =T(b u) + 3 u13 to be rewritten as 21 u + 2 - 1|
Introducing the scaled dual variable m = 2, the scaled form of the augmented Lagrangian becomes
L, u. m) = (bb + g(u) + -lbb- u + m
2 1mmll
(20)
This scaled form is better suited for implementing ADMM and SCPRSM schemes with proximal opera- tors (Parikh & Boyd 2013). The method of multipliers for 20 can be written
(k+1)) (k+1)) = argmin (b u (1) m (k)),
(21)
m (k+1) = m (k) (k++1) u (k+1).
(22)
The method of multipliers is generally not an implementable method since the primal update step 21 can be as hard to solve as the original problem [Bekk] 2017 Boyd et al.] 2011]. To overcome this challenge, ADMM employs an itcrative approach in the primal update step. In this approach, b and u are updatod sequentially in an alternating fashion, which is why the method called the alternating direction method of multipliers.
8

---

Published in Transactions on Machine Learning Research (11/2024)
An iterative scheme for the ADMM associated with 20 becomes
(k+1) argmin (b(b), (k) u (k) m 1
(23a)
b
u (k+1) := argmin L, (btt), (k) m (
(23b)
m (k+2) = m (k) + b (k+t1) (k+
(23c)
5.2 SCPRSM Framework
The difference between ADMM and PRSM in terms of convergence can be explained through the contraction properties their iterative sequences. The iterative sequence generated by ADMM is strictly contractive with respect to given solution set, whereas the sequence generated by PRSM contractive, but not strictly contractive (ee al. in PRSM, He et al. (2014) proposed incorporating relaxation factor > 0 into the Lagrange multiplier update steps, thus developing stricily contractive Peaceman--Rchhorre splitting method (SCPRSM). This modification ensures that the iterative sequence becomes strictly contractive, improving convergence prop- erties. The studics show that SCPRSM outperforms that ADMM generally leads to faster convergence
2014 Corman & Yuan 2014 He et al. 2002). To address the lack of strict contraction
compared to ADMM (Li & Yuan 2015
et al. 2021). The iterative scheme of the SCPRSM associated
with the augmented Lagrangian function 20
written as
(k-11) = argmin T, (b) (k), m m()))
(24a)
m (++ m (R) rbb(k+))
(k+1) = (k)),
(24b)
u (k-1) = argmin T,, (b (k+]), (k) m (k ())
(24c)
m (k-1) = m (++ + (b((++1) u (k+1)),
(24d)
where the parameter E (0,1) is relaxation factor. In addition to relaxation factor т, important distinction between SCPRSM and ADMM is the presence of an intermediate update for the multipliers, represented as m (*++) in the SCPRSM scheme. This step ensures balanced handling of the vectors b and
u, thereby leading to contractive iteration sequence that guarantees convergence to the solution of the
original optimization problem [H et al. 2014 Li Yuan 2015 Peaceman & Rachford, 1955). The iterative 23ch can be further simplified as
scheme of ADMM (23
(k++1) = argmin (f(b(() + 2
? 1lb(b)) (x) m (1))
u (k+1) argmin (g(u ()) + [(e(() - -(b(++1) (1)))
(25)
m (k+1) m + (++1) (1+
With proximal operators, we can now rewrite 25 as
(k++1) prox (n m (k)),
u (k+1) = prox, b (b+1) + m7 (k)),
(26)
m (k+1) = (k) b (x+1) 2 (l+11)
Similarly, the proximal version of Equations 24a
24d can be written
(k+1) = prpxooxxx yr (k) m (t));
m (R+ + = =m()) (bb(++1) u (k)). u (k+1) = proxx (b(++)) + m (++
(27)
m (k+1) = (+++2) + + (++1) - (x++)).

---

Published in Transactions on Machine Learning Research (11/2024)
5.3 Convergence Analysis
Now, we delve into the convergence properties of the iterative schemns of ADMM (23 and SCPRSM (24) This analysis will elucidate the conditions under which these methods converge and the nature of the solu- tions they yield. Specifically, the convergence of the proposed method is cstablished through the following theorems.
Theorem Let the sequences (b(h), (u( m ( be generated by the ADMM scheme 23c
23c
and its
Lagrangian is given by 20 Then the following three conditions hold:
(a) Sufficient decrease condition: For each iteration step k, _d1 > such that
(b(b(t))), (k+1), ml++l)) Lb(bb)) (k) m (k)) - 6(00
(b) Boundness condition: The sequences {b(h), (k), m (8) 10 are bounded and its Lagrangian (b(b),, (k) ml))) is lower bounded.
(c) Convergence: The Lagrangian in (20 is sponding sequence {b(M), (()), m²) converges to unique stationary point (b(*) u((), ml)))
(KL) function, then the corre-
Note that the function f(b) 11 Xbl is continuously differentiable function with respect to b. Its Then the Lipschitz constant for the gradient of the
gradient is computed as ff(b) = -x X6.
function f(b) can be computed as
V ( ((((+)))) s x -
< 2Amax (X x) b(++111
- 1(++
(28)
where || is the largest eigenvalue of x"X computed as Amax (X (X) We denote the Lipschitz X). gradient constant as l = 2Amaxx (x)X). The partial derivative of the Lagrangian 20 with respoct to b is given by
OgL,, u. m) (bbb + 7( u m)).
(29)
0 LL (b, m) 2x X + which is positive definite. This implics that 062
and the second partial derivative is the Lagrangian is strongly convex with respect to b. We will frequently use the propertics of the Lipschitz gradient constant of f(b) and and the strong convexity of (b, u, m) with respect to b in the proof (see Appendix B.2
Thcorem 2 Let the sequences (b(k), (k), m (k) be generated by the scheme grangian is given by 20 Then conditions (a)-(c) in Theorem ! hold.
24a
24d and its La-
The proof of Theorem 2 is found in Appendix
B.3
In conclusion, the sufficient decreasing and boundedness conditions are satisfied when the learning rate max{ 213 ./) in both Theorems 1 and 2 In practice, choosing involves a trade-off that requires 22 -/
careful consideration (see Section 5.5 The sufficient decreasing condition can verified for the mean squared error (MSE) loss because the update step inherently minimizes the loss, ensuring it decreases as the number of iterations increases.
5.4 Implementation
The key idea in the WL1LO optimization problem is split the variable = c d into two components: c, subject to L -norm regularization, and d, subject to L norm regularization. This decouples the optimization of and d within the loss function f(e d) |ll X(c + d)ll² allowing separate updates for each part.
10

---

Published in Transactions on Machine Learning Research (11/2024)
The docoupling is achieved through defining two translation functions, which manage the updates efficiently. The translation functions are defined as (u) = s(u + v) - a and Tu(()) = s(v + u) u. These enable alternating updates between C and d, using the current estimates of the other variable. By leveraging these translations, the loss function f(c + d) is effectively split, allowing the proximal operators to handle both I and - regularizations. Hence, for WLMLO-ADMM, the updates are made six steps, alternating between the two primal variables and v, with corresponding dual variables m and w. The steps are:
c(k+1) = prox (z) m (4),
u (A+1) prax, prax,y(c (c + +m (+)) (k)
m (k+1) = (A) + cli+1) - 1ktt)
(30)
d(++1) proxT
w (t))))
o ((+1) 1) = prox (d (1( + (8))) (k+1) (x) + ( (k+
For WLILO-SCPRSM, the updates made in eight steps
(k+1))= proXT, (u (z) m 1k)),
m 2) (b))
(b) + ((c - u (A)),
u (k+1) = prox,9 (c++)) + m (+++)),
m (k+1) 1) =m (l+ 1) ((c((t+11) r (kt.)),
(31)
(v)6 (o(x) - w (x)), w (K+ 1 = (A) + (+d()++)) - -v (A)))
(k+ 1) = pr)x
v (k+1) pr (d (k+) + (x(++1)),
w (k+1) == ( 1) + r (dd (k+12)).
Here, prox (c+ + mm) is the proximal operator for L', which is the soft- soft-thresholding function with learning rate y defined as
prox gy (c + m) (, + m) = max(0, |c + ml - 7)gg((c + m),
(32)
and pr + w) H (d + w) is the proximal operator for L which is hard thresholding operator defined as
d + w < v25,
H Hz6 (d + w)
d + w, if d w| > v25, {0, d + w
(33)
if |d + w| = /25.
The iterations are terminated when convergence is reached according to (k(c(t) + dd()) (z) + (l)) 5 3(1 + (lm + (+)))))>>) for tolerance parameter which was set to 10-5
For comparison purposes, also implement the LASSO, SCAD, MCP and EN methods using the proximal ADMM and SCPRSM schemes (see Appendix D Note that the convergence analysis of (30) and (31 is straightforward from 23 and 24 However, the introduction of translation functions and the variables e and d increases the dimensionality of the optimization problem, making the thcoretical analysis very extensive. Therefore, we omit the detailed proof of the algorithm that involves the translation functions.
5.5 Determining the Learning Rate
Choosing the learning rate (step size) is crucial for efficiency and proper convergence of optimization algo- rithms. There are two main methods for determining the learning rates y and 8 (Bock 2017 Bertsekas 2016 Boyd & Vandenberghe, 2004): 1. Backtracking line-scarch: This method adjusts the learning rate iteratively based specific criteria. However, it is both computationally expensive and time-consuming, it requires
11

---

Published in Transactions on Machine Learning Research (11/2024)
multiple evaluations of the objective function and its gradient during the search process. These repeated evaluations increase the overall computational load, particularly in high-dimensional problems. 2. Constant learning rate: In contrast, this method uses a fixed learning rate throughout the whole optimization process. It is computationally simpler and avoids the time overhead associated with frequent adjustments, making it more efficient in many scenarios.
We adopt the constant learning rate approach, using 7 1 la on the training set defined as
for all K, where ||||| is the operator norm
|X||2 max ||b1
||
(34)
Equivalently, || 22 is the maximum singular value of X (amax (X)) which measures the maximum amount by which the matrix X can stretch a vector relative to its original length (Horn & Johnson 2013). The supremum-based definition can be applied in (34), particularly in infnittediimeniional contexts (Bhatia
1997). However, in finite dimensions, the maximum and supremum coincide for the operator norm defined in (3-h. By setting the learning rale to 2 we ensure that the learning rate is scaled appropriately relative
to the maximum possible stretch of X. We use the same formula for 8.
5.6 Bayesian Optimization for Hyperparameter Tuning
Tuning the regularization parameter X, the weight parameter a and the relaxation factor r via cross-validation or grid search can be computationally expensive. Bayesian Optimization (BO) is a more advanced, data- driven approach which offers probabilistic model-based method for hyperparameter tuning (Gao et al.
2021 Shahriari ct al., 2015). For the latest advancements, see Wang al. (0223 and Yang ct al. (2024).
BO uses surrogate model, often a Gaussian Processes (GP), to approximate the true objective fune- tion. Hyperparameters are collected in & = [a, A. and the objective function <[d] is modeled as d0] ~ GP(m(U), A(0, I')), where m|U] is its mcan and k(U, 0'] the kernel (variance) function. The objec- tive function is evaluated at j sequential points MSE = L(06)),, with MSE ~ N(((00)), o2)) This process induces a posterior over the acquisition function, guiding the selection of the next hyperparameters. Common acquisition functions include probability of improvement (PI), expected improvement (EI), upper confidence bound (UCB), and mutual information (MT) (Snoek et al., 2012). BO starts with an initial set of hyperparameters and objective function values to train the surrogate model. The acquisition function balances the postcrior mean (1)) for exploitation and variance (v(1) for exploration. The GP-UCB is given by
Φ0-1) aargmax(=(0)) + xv())),
0
where (1) is driven by the mean function m(),, v(0) by the variance function k(U), and x de- termines the trade-off between exploitation and exploration. Contal et al. (204)) improved GP-UCB with the Gaussian Process Mutual Information algorithm (GP-MT) as 00-1)
(gmaa((((0)) +
vlog(2/e)( E(16)) +5(1-1)- Vs(-Ui)), where ş controls exploration, < e < 1, and E(1() is the variance function at 861)
6
Numerical Experiments
6.1
Materials
We evaluate our proposed method using simulated genomic dataset well two real-world genomic datasets. Specifically, meclotie polymorphisms (SNPs), which is type of genetic variation that represent differences in one of the two muclcotides that make up an individual's DNA at a specific location compared to the most common nucleotide pair found in a population. SNPs are typically represented by a count of 0, 1, or 2, where means both nucleotides at the SNP location match the most common pair, 1 indicates that one of the two nuclootides differs from the common pair and 2 mcans both nuclcotides at the location differ from the most common pair. This count system helps quantify genetic variation at a specific SNP site. A detailed explanation of these datasets is provided in Appendix C
12

---

Published in Transactions on Machine Learning Research (11/2024)
Simulated QTLMAS 2010 Dataset (Szydlowski & Paczyúskka 2011 p: This dataset comprises 3226 individuals,
with genomic single nucleotide polymorphism (SNP) data organized in matrix X of size 3226 x 9723, with two observed traits (response variables): quantitative and a binary trait. In this study, the quantitative trait was chosen as the phenotype which is represented as a continuous response vector (y) of length 3226.
Real Pig Dataset (Cleveland et al. 2012 : This dataset contains genomic SNP data from 3534 individuals, organized in a matrix of size 3534 52842, along with phenotypic data for five traits. We used trait 4, which had a heritability of 0.58, as the phenotype that constitutes a vector (y) of length 3534.
Real Mice Dataset (Pérez & de Los Campos 2014): This dataset contains data from 1814 individuals, with a genomic SNP data matrix of size 1814 x 10346, along with two traits: body length (BL) and body mass index (BMI). In this study, the continuous trait BL was chosen to be our response vector (y) that has a length of 1814.
6.2 Results
The WLIL0-ADDM, WLILO-SCPRSM, EN-ADMM, EN-SCPRSM, LASSO-ADMM and LASSO-SCPRSM methods were implemented in Julia 1.10.1 (Bezanson et al. 2017) using the ProximalOperators package (An- tonello al., 2018). For SCAD-ADMM, SCAD-SCPRSM, MCP-ADMM and MCP-SCPRSM, we wrote our own code manually in Julia. For all methods, the BO was performed with the BayesianOptimization package using an ElasticGPE model and the squared exponential automatic relevance determination (SEArd) ker- nel (Firbrotther al. 2018). The initial values of b, è and d were set to the marginal covariances bet ween y and X, multiplied by 0.0001. By conducting preliminary runs for each set of hyperparameters using BO, we identified the optimal range of parameters. BO with the MI acquisition function was executed for hy- perparameter tuning of all methods. The regression coefficients of the model are obtained from the training dataset, and once the model is trained, it predicts outcomes on the test dataset. The MSF is then calculated on the test dataset to assess the model's generalization performance. The test MSE was monitored during the BO process to ensure convergence, which was indicated by no further decrease in MSE. All analyses were executed on a Linux computing platform equipped with an AMD EPYC 7302P 16- 16-Core Processor and 32GB of system memory.
6.2.1 Simulated QTLMAS 2010 Dataset
BO was executed for 250 iterations with GP function evaluations per iteration across all methods. The lower and upper bounds for X1 were set to 0.001 and 1000.0, 0.1 and 1800.0, and 0.001 and 600.0 for LASSO- ADMM, SCAD-ADMM, and MCP-ADMM, respectivcly. The lower and upper bounds for r were set to 0.01 and 0.999, 0.001 and 1.0, 0.001 and 1.0, and for λι were set to 0.001 and 1000.0, 0.1 and 1800.0, and 0.001 and 600.0 for LASSO-SCPRSM, the lower and upper bounds for λ were set to 10.0 and 600.0 and for X2, they were set to 0.001 and 1.0, respectively. For EN-SCPRSM the lower and upper bounds for X were set to 10.0 and 500.0, for X2 they were set. 0.001 and 200.0, for r, they were set to 0.01 and 0.99 respectively. For WL.M the lower and upper bounds for a were set to 0.01 and 0.99, and for ' they were set to 0.001 and 500.0, respectively. For WLILO-SCPRSM, the lower and upper bounds for were set to 0.0001 and 0.999, for r they were set to 0.0001 and 1.0, and for Xı they were set to 0.0001 and 500.0, respectively. The best result, with minimum test MSE of 64.55, was found with W.ILO-SCPRSM at X = 1991.5,, = 0.90, and r = 0.48 (Table I. Timing of the last evaluation with optimized parameters showed that SCAD-ADMM executed most quickly in only 6.68 seconds. It should be noted that those methods with one regularization parameter tend to be faster to train compared to other methods with two or three hyperparameters.
and MCP-SCPRSM, respectively. For EN-ADMM
6.2.2 Real Pig Dataset
For the Pig dataset, we employed 5-fold cross-validation with random allocations into training and test data to obtain the minimum test MSE on the test data set, with the results averaged over the folds. Here, for all methods, BO was executed for 100 iterations with GP function evaluations per iteration due to the large dataset size. The lower and upper bounds for X were set 50.0 and 600.0, 0.1 and 1000.0, and 0.1 and 20.0 for LASSO-ADMM, SCAD-ADMM, and MCM-ADMM, respectively. The lower and upper bounds for
13

---

Published in Transactions on Machine Learning Research (11/2024)
Method
min MSE
À1
X2
Time (s)
Number of non-zeros
LASSO-ADMM 66.55
294.22
9.31
417
LASSO-SCPRSM
65.95
312.51
0.19
21.51
334
SCAD-ADMM
66.50
309.65
6.68
386
SCAD-SCPRSM
65.91
299.79
0.42
19.68
352
MCP-ADMI 69.89
362.57
10.10
3898
MCP-SCPRSM
68.12
268.17
0.99
17.93
3883
EN-ADMM
66.52
307.65
0.79
11.87
390
EN-SCPRSM 65.92
288.15
0.001
0.99
22.18
375
WLILO-ADMA 64.77
370.46
0.86
28.63
324
WLILO-SCPRSM 64.55
391.55
0.90
0.48
25.39
275
Table 1: Performance evaluation of various regularization methods with optimal parameters on simulated QTLMAS data. The best-performing test MSE and most sparse model are highlighted in bold.
A were set to 50.0 and 400.0, 50.0 and 400.0, and 0.01 and 20.0, and lower and upper bounds for r were set to 0.01 and 1.0, 0.01 and 1.0, and 0.01 and 1.0 for LASSO-SCPRSM, SCAD-SCPRSM, and MCP-SCPRSM, respectively.
For EN-ADMM the lower and upper bounds for were set 10.0 and 600.0 and for d2, they were set to 0.001 and 1.0 respectively. For EN-SCPRSM the lower and upper bounds for ê were set to 0.1 and 200.0, for X2, they were set to 0.01 and 100.0, for r, they were set to 0.001 and 1.0, respectively. For WL.1L0-- ADMM the lower and upper bounds for were set to 0.001 and 0.99, and for X1, they were set to 0.001 and 100.0, respectively. For WLILO-SCPRSM, the lower and upper bounds for a were set to 0.001 and 0.99, for they were set to 0.001 and 1.0, and for A they were set to 0.0001 and 200.0, respectively, We obscrved little variability the minimum test MSE across the CV-folds for all methods. Hence, we report the mean minimum test MSE using the average estimates of the respective parameters for all methods. The best result, with mean minimum test. MSF of 4.48, found with WT.1.0-SCPSS. with mean estimates A 240.63, 0.54, and 0.41 (Table 2). The average timing over the folds of the last evaluation with was fastest, taking only 25.8 seconds.
optimized regularization parameters showed that
Method min MSE
A1
A2
Time (s)
Number of non-zeros
LASSO-ADMM
4.53
118.75
26.2
1515
LASSO-SCPRSM
4.50
115.63
0.32
48.21
1200
SCAD-ADMM
4.64
127.80
25.8
1272
SCAD-SCPRSM
4.50
115.63
0.32
48.68
1200
MCP-ADMM
6.13
100.0
27.87
22909
ACP-SCPRSM
6.12
106.25
0.94
25.58
21627
EN-ADMM
4.53
120.63
0.31
26.86
1447
EN-SCPRSM
4.51
118.79
96.88
0.97
32.42
1433
WL.11.-ADMM
4.49
43.75
0.56
130.82
1093
WLILO-SCPRSM 4.48
240.63
0.54
0.41
124.89
852
Table 2: Evaluation of the performance of various regularization methods with optimal parameters, averaged across five CV-folds on the pig data. The best-performing Lest. MSE and most sparse model are highlighted in bold.
6.2.3 Real Mice Dataset
Similar to the Pig datasel, we employed 5-fold cross-validation also for this data. BO was executed for 100 iterations with GP function evaluations per iteration across all methods. The lower and upper bounds for X, were set 0.0001 and 20.0, 0.001 and 20.0, and 0.1 and 20.0 for LASSO-ADMM, SCAD-ADMM, and respectively. The lower and upper bounds for A1 were set 0.001 and 35.0, 0.001 and 35.0, and 0.01 and 20.0, and lower and upper bounds for set to 0.001 and 1.0, 0.001 and 1.0, and 0.01 and
14

---

Published in Transactions on Machine Learning Research (11/2024)
1.0 for LASSO-SCPRSM, SCAD-SCPRSM, and MCP-SCPRSM, respectively. For EN-ADMM the lower and upper bounds for Å were set to 0.01 and 18.0 and for X2, they were set to 0.001 and 1.0 respectively. For EN-SCPRSM the lower and upper bounds for Y were set to 0.1 and 42.0, for 22, they were set to 0.01 and 40.0, for r, they were set to 0.01 and 1.0, respectively.
For WLIL.-ADMM, the lower and upper bounds for a were set 0.001 and 0.99, and for Å1, they were set to 0.001 and 100.0, respectively. For WLILO-SCPRSM, the lower and upper bounds for were to 0.001 and 0.99, for they were set to 0.001 and 1.0, and for X they were set to 0.0001 and 200.0, respectively. The best result, with mean minimum test MSE of 0.259 was found with WLLLO-SCPRSM at the average estimates \ 56.25, a = 0.28, and r= 0.16 (Table 3).
The average liming over the folds of the last evaluation with optimized regularization parameters showed that SCAD-ADMM once again was fastest with a time of only 2.08 seconds.
Method LASSO-ADMM
min MSE
A1
A1
Q
r
Time (s)
Number of non-zeros
0.273
20.0
2.10
319
LASSO-SCPRSM
0.267
24.06
0.81
3.14
280
SCAD-ADMM
0.274
24.0
2.08
319
SCAD-SCPRSM
0.267
24.06
0.81
2.78
280
MCP-ADMM
0,273
20.0
2.14
432
MCP-SCPRSM
0.273
20.0
0.01
26.62
387
EN-ADMM
0.276
18.0
0.99
2.15
539
EN-SCPRSM 0.267
24.98
38.75
0.97
2.51
256
WLILO-ADMM 0.265
43.75
0.56
19.97
216
WLILO-SCPRSM 0.259
56.25
0.28
0.16
22.83
188
Table 3: Evaluation of the performance of various regularization methods with optimal parameters, aver- aged across five CV-folds on the mice dataset. The best-performing test MSE and most sparse model are highlighted in bold.
7 Discussion
The WL1LO method demonstrates superior performance across all datasets by achieving the lowest. MSE and the fewest non-zero coefficients. This highlights its effectiveness and efficiency as regularization technique in high-dimensional data analysis, making it a valuable altornative to the LASSO, SCAD, MCP and EN. The weighting parameter a in WL1LO provides flexibility tuning the regularization effect, making the method adaptable to different datasets and problem settings. This adaptability enhances robustness and applicability across diverse scenarios.
It has been demonstrated several times that the SCAD and MCP often outperform the LASSO (aan et al, 2014a; Fan & Li 2001 Zhang 2010). However, while the LASSO, SCAD, MCP and EN also offer competitive approaches to regularization, the WLILO method consistently outperforms them, providing enhanced model sparsity and interpretability without compromising predictive accuracy. The joint sparsity induced by the I, and the 7.0 components make the resulting model more interpretable. This crucial many scientific and industrial applications, where understanding the model is important as its prodictive power.
The use of the SCPRSM algorithm introduces additional parameter r. which allows for finer control over the optimization process and potentially leads to better convergence properties and more precise model fitting. Across all datasets used, the SCPRSM variants demonstrate strong performance by achicving the smallest minimum MSEs while maintaining a manageable number of non-zero coefficients. Specifically, WLILO-SCPRSM consistently achieves the lowest MSE across all our datasets, demonstrating its superior ability minimize prediction errors. is likely that will be highly effective terms of both accuracy and reliability across other types of data. Several other studies have shown that SCPRSM outperforms ADMM [Li & Yuan 2015
Li et. al. 2021
15

---

Published in Transactions on Machine Learning Research (11/2024)
In this paper, we have mostly focused the regularization part and note that there certainly is room for computational advancements. ADMM can be improved using techniques such accelerated ADMM (Zhang) et al. 2019 Zeng et al. 2024 and references therein) as well as stochastic distributed ADMM (Chen et al. 2021 and references therein). Similarly, for SCPRSM, further computational improvements can be achievod
via stochastic SCPRSM (Na el al. 2017) and indefinite-proximall SCPRSM (Gu al. 2022 Bai el al. 2023
8 Conclusion
This paper introduces a novel joint weighted L'- and L' method denoted WL1LO based on proximal mappings and translation functions, aiming to debias the bias introduced by the L'-morm when applied to high-dinensional data. Our model introduces weighting parameter a, allowing for the adjustment of the influence of both regularizers. The convergence of ADMM and SCPRSM for the developed method is shown under reasonable assumptions. All hyper are optimized using Bayesian optimization. The WLILO-SCPRSM mothod consistently achicves the lowest MSE across all datascts when compared to all other tested regularization methods (LASSO, EN, SCAD and MCP). Hence, the WLILO-SCPRSM's superior performance across different genomic high-dimensional datasets demonstrates its versatility. Our current paper focuses primarily on prediction. In future work, we plan to specifically address the properties of variable selection.
Acknowledgements
We acknowledge funding from the University of Oulu & the Academy of Finland Profi 326291.
References
Niccolò Autonello, Lorenzo Stella, Panagiotis Patrinos, and Toon Van Waterschoot. Proximal gradient algorithms: Applications in signal processing. arXiv preprint arXiv: 180300000.11 2018.
Hédy Attouch, Jérôme Bolte, Patrick Redont, and Antoine Soubeyran. Proximal alternating minimiza- tion and projection methods for nonconvex problems: An approach based on the incquality. Mathematics of Operations Research, 35(2):438-457, 2010.
Hedy Attouch, Jérôme Bolte, and Benar Fux Svaiter. Convergence descent methods for semi-algebraic and tame problems: Proximal algorithms, forward backward splitting, and regularizod Causs-Seider methods. Mathematical Programming, Series A, 137(1):91 129, 2013.
Jian-Chao Bai, Feng-Miar Bian, Xiao-Kai Chang, and Lin Du. Accelerated stochastic peaceman-rachfon method for empirical risk mininization. Journal of the Operations Research Society of China, 11(4): 887-807,40 2023.
Amir Beck. First-order Methods in Optimization. SIAM, 2017.
Dimitri Bertsekas. Nonlinear Programming, 3rd ed. Athena Scientific, Nashua, NH, 2016.
Dimitris Bertsimas, Jean Pauphilet, and Bart Van Parys. Sparse regression: Scalable algorithms and empir- ical performance. Statistical Science, 35(4):555 578, 2020.
Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. Julia: A fresh approach numerical computing. SIAM Review, 59(1):65 98, 2017.
Rajendra Bhatia. Matriz Analysis, volume 169. Springer-Vett New York, 1997.
Jérôme Bolte, Shoham Sabach, and Mare Teboulle. Proximal alternating linearized minimization for non- convex and nonsmooth problems. Mathematical Programming, 146(1):459944)) 044449994444444 2014.
Stephen Boyd, Neal Parikh, Eric Chu, Borja Pelcato, Jonathan Eckstein, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends@ in Machine Learning, 3(1):1-122, 2011.
16

---

Published in Transactions on Machine Learning Research (11/2024)
Stephen P Boyd and Lieven Vandenberghe. Conver Optimization. Cambridge University Press, 2004.
Peter Bühlmann and Sara Van De Geer. Statistics for High-Dimensional Data: Methods, Theory and Applications. Springer-Verlag Berlin Heidelberg, Berlin, Germany, 2011.
Hao Chen, Yu Ye, Ming Xiao, Mikael Skoghund, and H Vincent Poor. Coded stochastic admm for decen-
tralized consensus optimization with edge computing. JEEE Internet of Things Journal, 8(7):5300-5373, 2021.
Victor Chernozhukov, Christian Hansen, and Yuan Liao. lava attack on the recovery of sums of dense and sparse signals. The Annals of Slatistics, 45(1):39 76, 2017.
Francis H Clarke, Yuri S Ledyaev, Ronald J Stern, and Peter R. Wolenski. Nonsmooth Analysis and Control Theory, volume 178. Springer Science & Business Media, 2008.
Matthew A Cleveland, John M Hickey, and Selma Forni. A common dataset for genomic analysis of livestock populations. G3: Genes/ Genomes! Genetics, 2(4):429-4355 2012.
Emile Contal, Vianney Perchet, and Nicolas Vayatis. Gaussian process optimization with mutual information. International Conference on Machine Learning, volume 32, pp. 253 261. PMLR, 2014.
Etienne Corman and Xiaoming Yuan. A generalized proximal point algorithm and its convergence rate. SIAM Journal on Oplamization, 24(4):1614 1638, 2014.
Jamie Fairbrother, Christopher Nemeth, Maxime Rischard, Johanni Brea, and Thomas Pinder. Gaussian- processes. jl: A nonparametric Bayes package for the Julia language. arXiv preprint or 181 2.209904. 2018.
Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the Americon Statistical Association, 96(456):1348 1360, 2001.
Jianqing Fan and Jinchi Lv. A selective overview of variable selection high dimensional fcature space. Slatistica Sinica, 20(1):101-11) 148, 2010.
Jianqing Fan, Jinchi Lv, and Lei Qi. Sparse high dimensional models in economics. Annual Review of Economics, 3(1):291-311 -317, 2011.
Jianqing Fan, Yingying Fan, and Emre Barut. Adaptive robust variable selection. Annols of Slatistics, 42 (1):324-35)) 2014a.
Jianqing Fan, Fang Han, and Han Lin. Challenges of big data analysis. National Science Review, 1(2): 293 314, 2014b.
Masao Fukushima and Hisashi Mine. A generalized proximal point algorithm for certain non-convex mini- mization problems. International Journal of Systems Science, 12(8):989 1000, 1981.
Haiping Gao, Shifa Zhong, Wenlong Zhang, Thomas Igou, Eli Berger, Elliot Reid, Yangying Zhao, Dylan Tambeth, Lan Gan, Moyosore A. Afolabi, et al. Revolutionizing membrane design using machine learning- Bayesian optimization. Encironmental Science & Technology, 56(4):2572 2581, 2021.
Christophe Girand. Introduction to Highi-Dimensions Statistics. CRC Press, Boca Raton, FL, 2015.
Yan Gu, Bo Jiang, and DR Han. An indefinite- petttte--rooiim---aaeeed strictly contractive peaceman-racchord split- ting method. Journal of Computational Mathematics, 41:1017 1040, 2022.
Trevor Hastic, Robert Tibshirani, and Martin Wainwright. Statistical Learning with Sparsity: The Lasso and Generalizations. CRC Press, New York, 2015.
Trevor Hastie, Robert Tibshirani, and Ryan Tibshirani. Best subset, forward stepwise or Lasso? Analysis and recommendations based extensive comparisons. Statistical Science, 35(4):579 -122, 2020.
17

---

Published in Transactions on Machine Learning Research (11/2024)
Bingsheng Hc, Li-Zhi Liao, Deren Han, and Hai Yang. A new inexact alternating directions method for monotone variational inequalities. Mathematical Programming, 92:103-118, 2002.
Bingsheng He, Han Liu, Zhaoran Wang, and Xiaoming Yuan. A strictly contractive Peaceman Rachford splitting method for convex programming. STAM Journal Optimization, 24(3):1011-1044) 2014.
Georg Heinze, Christine Wallisch, and Daniela Dunkler. Variable selection a review and recommendations for the practicing statistician. Biometrical Journal, 60(3):431-499, 2018.
R. R. Hocking and R N Leslie. Selection of the best subset in regression analysis. Technometrics, 9:531-544 1967.
Arthur E Hoerl and Robert W Kennard. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1):55-67, 1970.
Roger Horn and Charles R Johnson. Matrix Analysis. Cambridge University Press, 2013.
Tain M Johnstone and Michael Titterington. Statistical challenges of high-dimensional data. Philosophical Transactions of the Royal Society A: Mathematical, Physiiaal and Engineering Sciences, 4253, 2009.
Alexander Kaplan and Rainer Tichatschke. Proximal point methods and nonconvex optimization. Journal of Global Optimization, 13:389 106899406,,, 1998.
Hai Yen Le. Generalized subdifferentials of the rank function. Optimization Lellers, 7:731 743, 2013.
Peixuan Li, Yuan Shen, Suhong Jiang, Zchua Liu, and Caihua Chen. Convergence study on strictly contrac- rachford splitting method for nonseparable convex mininization models with quadratic coupling terms. Computational Optimization and Applications, 78:87- 124,124, 2021.
live peaceman rachford
Xinxin Li and Xiaoming Yuan. A proximal strictly contractive convex programming with applications imaging. SIAM Journal on Imaging Sciences, 8(2):1332 1365, 2015.
splitting method for
Xingran Liao, Xuekai Wei, and Mingliang Zhou. Minimax concave penalty regression for superresolution image reconstruction. IEEE Transactions on Consumer Electronics, (0(1):2999-300 90099 2023.
Christos Louizos, Max Welling, and Diederik P Kingma. Learning sparse neural networks through lo regu- larization. arXiv preprint 11t::1122.01111111 2017.
Chris J Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. arXiv preprint ar arXiv: 111111111001111 2016.
Boris S Mordukhovich. Varialional Analysis and Generalized Differentiation II: Applications, volume 331. Springer, 2006.
Sen Na, Mingyuan Ma, Shuming Ma, and Guangju Peng. Stochastic strictly contractive peaceman-rachforl splitting method. Xrxiv preprint 3XXx::7111.0455 2017.
Joseph O Ogutu and Mans-Peter Piepho. Regularized group regression methods for genomic prediction: Bridge, MCP, SCAD, group bridge, group lasso, sparse group lasso, group MCP and group SCAD. BMC proceedings, volume 8, pp. 1 9. Springer, 2014.
Neal Parikh and Stephen Boyd. Proximal algorithms. Foundations and Trends in Optimization, 1(3): 123 231, 2013.
Donald W Peaceman and Henry H Rachford, Jr. The numerical solution of parabolic and clliptic differential equations. Journal of the Sociely for Industrial and Applied Muthematics, 3(1):28-44 41, 1955.
Paulino Pérez and Gustavo de Los Campos. Genome-wide regression and prediction with the BGLR statis- tical package. Genetics, 198(2):483-444 -45,, 2014.
18

---

Published in Transactions on Machine Learning Research (11/2024)
Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando De Freitas. Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, 104(1):148-17,, 2015.
Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical Bayesian optimization of machine learning algorithms. Aduances in Neural Information Processing Systems, 25, 2012.
Macicj Szydlowski and Paulina Paczyńska. QTLMAS 2010: Simulated dataset. In BMC procecdings, vol- ume 5, pp. 1 3. Springer, 2011.
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Melhodological), 58(1):267 288, 1996.
Laura Tologi and Thomas Lengauer. Classification with correlated features: unreliability of feature ranking and solutions. Bioinformaliees (7114):1986-1999 2011.
Martin J Wainwright. High-dimensional Statislics: A Asy-poopie Vieupoinl. Cambridge University
Press, Cambridge, United Kingdom, 2019.
Patrik Waldmann. A proximal LAVA method for genome-wide association and prediction of traits with mixed inheritance patterns. BMC Bioinformatics, 22(1):1 16, 2021.
Hao Wang, Zhanglei Shi, Chi-Sing Leung, and Hing Cheung So. ADMM-MCP framework for sparse recovery with global convergence. IEEE Transactions on Signal Processing, 2018.
Ting Wang and Hongwei Liu. class of modified accelerated proximal gradient methods for nonsmooth and nonconvex minimization problems. Numerical Algorithmss, 95(1):207-22) 241, 2024.
Xilu Wang, Yaochu Jin, Sebastian Schmitt, and Markus Olhofer. Recent advances in Bayesian optimization.
ACM Computing Surveys, 55(13)):1-30 2023.
Yu Wang, Wotao Yin, and Jinshan Zeng- Global convergence of ADMM in nonconvex nonsmooth optimiza- tion. Journal of Scientific Computing, 78:29-63, 2019.
Yuqia Wu, Shaohua Pan, and Shujun Bi. Kurdyka-tojasioviez property of zero-norm composite functions.
Journal of Oplimization Theory and Applicutions, 188:94 112, 2021.
Kaixin Yang, Long Liu, and Yalu Wen. The impact of Bayesian optimization on feature selection. Scientific Reports, 14(1):3948, 2024.
Maryam Yashtini. Convergence and rate analysis of a proximal linearized admm for nonconvex nonsmooth optimization. Journal of Global Oplimization, 84(4):913-99) 939, 2022.
Jihun Yun, Peng Zheng, Eunho Yang, Aurelie Lozano, and Aleksandr Aravkin. Trimming the & regular-
izer: Statistical analysis, optimization, and applications to deep learning. In International Conference Machine Learning, pp. 7242 7251. PMLR, 2019.
Jihun Yun, Aurélic C Lozano, and Eunho Yang. Adaptive proximal gradient methods for structured neural networks. Advances in Neural Information Processing Systems, 34:24365-24338 2021.
Yuxuan Zeng, Zhiguo Wang, Jianchao Bai, and Xiaojing Shen. An accelerated stochastic stochassti and nonsmooth finite-sum optimization. Automatica, 163:111554, 2024.
admm for nonconvex
Cun-Hui Zhang. Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics, 38(2):894 942, 2010.
Juyong Zhang, Yue Peng, Wenqing Ouyang, and Bailin Deng. Accelcrating admm for efficient simulation and optimization. ACM Transactions on Graphics (TOG), 38(6): 21, 2019.
Peng Zhao and Bin Yu. On model selection consistency of Lasso. The Journal of Machine Learning Research, 7:2541-22 1-663, 2006.
19

---

Published in Transactions on Machine Learning Research (11/2024)
Tuo Zhao, Han Liu, and Tong Zhang. Pathwise coordinate optimization for sparse learning: Algorithm and theory. The Annals of Slatistics, 46(1):180-21) 2181800-211-000 2018.
Liu Ziyin and Zihao Wang. spred: Solving L1 penalty with SGD. In International Conference Machine Learning, pp. 43407-43422, PMLR, 2023.
Hui Zou and Trevor Hastie. Regularization and variable selection via the elastic net. Journal of the Royal Stalistical Sociely: Series B (Stalistical Melhodology), 67(2):301-33) -2000, 2005.
Appendix
A Related Problems
Consider the following formulation
b = argmin F(b) = [bb) + A(1 - a))l
(35)
where I(6) = 11y xb + Aa Aa|bbll Since r(bb) is a convex function, has a global minimun value. By the Weierstrass theorem, continuous function over nonempty compact set attains minimum. The existence of an optimal solution is guaranteed if a function is continuous over closed set and coercive over the set (Bertsekas] 2016). Beck (2017) demonstrates that the latter extends to closed functions, i.e. a closed and coercive function over a closed set attains an optimal solution.
For the case |16||o
Σ 1(bs + 0) with X(1 a) > 0, we need to show it is closed function. Let
g(b)= Σ
I(bi,, where 1: + 10.1} is defined
1=1
[A(1 a))-a),
+ 0, b. = 0.
r(b:)
0.
The function I(.) is closed since its level sets, given by
0,
1) < 0,
Lev(I,n)
{0},
17 [0. 1),
(36)
R.
2 1,
are closed sets. Here, g is a closed function. Furthermore, using Theorem 2.6 in (Beck 2017). the closedness of | implics its lower semi-continuity.
A vector b* is local minimum of the function F, if there exists 0 such that F(b") < F(b) for all b E R with b 1 b"ll < e. A vector 6* is a global minimum if F(b") < Fbb) for all b E IP.
For illustration purposes, we generated a random design matrix X with dimension 100 x 500 by simulating 100 samples and 500 features. For each value in the range between -5 and 5, the outcomes function F(b) that combines the squared error loss with the regularization term that consists of the weighted sum of the L and Zo-orsms were produced. The regularization parameters were set to À and 0.6. Therefore, the plotted F(b) includes both the error and regularization terms, rather than solely the penalty norms. One can see that F(b) is nonconvex because any point between the endpoints A and B, as indicated by the dashed red line in Figure I lies outside the domain of F(b). In fact, the shape of the function F(b) is similar
to that of nonconvex regularization methods such as SCAD and MCP (Fan & Li 2001 Zhang 2010 Zhao
et al. 2018
Various shapes of the function F(b) for different values of and À 1 are depicted in Figure to illustrate the impact of the function's behavior.
20

---

Published in Transactions on Machine Learning Research (11/2024)
3000
F(b)
1000
Figure 1: Ilustratio of the nonconvexity of the function F
2500
0.0
0.2
0.8
2000
1.0
1500
F(b)
1000
500
0
-2
b
Figure 2: Different function values of F' with regularizer A and 0 (22),02,,550088.111 (L').
B Proofs
B.1 Proof of Lemma 1
Given f(u) (m, u) + a, the proximal operator of the function f(u) is defined as prox,(u) argmin (f(v) + -1ll 1--ul wl12]. Substituting the affine function f(v) (m, v) + gives us prox,(u)
argmin (m, v) + a 1 u 1 The term a is constant and does not affect the minimization (it doesn't 't
depend on v). so we can ignore it. Next, the squared norm term Hll ulB can be expanded as el - u13 w 3 - - + [lul(). Since ||u| is constant with respect to v, it can also be ignored in the minimization. Now, combining the linear terms involving v, we have argmin ( 11 + (m u, }
This expression is quadratic function in To minimize it, we set the gradient of the objective function
21

---

Published in Transactions on Machine Learning Research (11/2024)
1 with respect to equal to zero. The gradient of =110112 is v. The gradient of (m u, v) is u. Set- 2 ting the total gradient zero gives us: v + (m u) O. Solving for v, find m. Therefore, pro u - m.
B.2 Proof of Theorem
While our model differs from that of Wang et al. 22019 2018), we adopt similar proof framework.
Proof (a): From [23a], ((++1) minimizes T,, (b,),, ( mlf) and since L, (b. u, m) is strongly convex with
respect to b, the Lagrangian function satisfies the following inequality [Beck 2017
(b((1)),, u (R) m (+)) (b(k), u (K) m (0) - - [1 (k+1)) 6 12.
(37)
2
where L-(-) is a p-strongly convex function (p 0). From the augmented Lagrangian function in (20). we have
(b)+)) (b(k+1,, (k+11, m ((t1)) - L, (b(k+1), u (1++1), -m m!! (
2 (blitt) - -u ulattl) (m - m ml())
(38)
Now we rewrite 23c as
blk+1) - ulk+1) m (k+1) mk
(39)
From
23b and 29
we get
(((b(t))) + >b((+t)) - u(k+1) + (b)) = O.
(40)
Substituting 23c into 40
we obtain
= =
(41)
Using 39 and 41
38
becomes
- 7 (l m 101) =7 ( F(o 1 F
17 (1b(kt))- -bw0)22)
- ?
(42)
Hence,
13 (b(k+1)), (k++1), m (k+1)) L, (b((+1), (k+1), m (k)) s (16th+2)- 1(()
(43)
Here, the term 1 2 denotes a Lipschitz gradient of the function f(b). From 23b we have that
m (+)) (b(H)) u (k), ml))))) s L,, (b(++1), uu (1+1),
(44)
Finally, combining
37 0
43 and
44 we obtain the desired inequality as
L, (b re (k+1) m (+++))) (b(h), .re m (8))
= (b((tt) (k+1), m (k+1) - L, r (1+1) m (k)
+ (bbltt) u (k+1) m (k)) (6(K2))) (t) m(())
+ L (b+1) (M) m (k)) (b(1), ee( m) (k)
12
| 1
<
2)
S1 b (M++)) - 13
(45)
P 12 where 0 2 7
212 and r Hence, the sufficient decreasing condition is met.
Proof (b): We utilize the descent lemma to prove that T. (b(k), (k) .m (t) is lower bounded for any k.
22

---

Published in Transactions on Machine Learning Research (11/2024)
Lemma 2 (Descent lemmaj Let the function f belong to the class of continuously differentiable functions with constant ly Lápschitz continuous gradients. Then for any two points ((k) and (8)
= s f(b(t) + V - b(k) + 4 I1 11.
(46)
2
The proof of the descent lemma can be found As result of the Descent lemma, the sequence is lower bounded as
Beck 2017), see Lemma 5.7.
I,((bh), (1) m (k) =f(b) + g(u) + 2/1 bb + - m - m mm =f(b())) + g(u (81) mm' TT (b((tt) () -x) + (T/2)
?
u 12
+ 7 11
2 (() + g(((k))
| - m m
2 2
2 0o for 2 1z.
(47)
Hence, from
47
Ly(b(h, ( ()) is lower bounded.
As established in the proof (a), the sufficient descent property implies that Ly (b(k), (k), ()) is upper- bounded by (b", u m²). To prove that the sequence (b(A,, (*)), ml is bounded, we start by rewriting
45p as
[b(++]] - 11 - 81
(b(k), _k!)) - I, (b(it1), (l++1), m (b+))))
[k(++)) - 11 v 81 (b", u. uu. m (b wtt1 /_timm 1)
(48)
Equation 48 also holds as l + 0o. Hence, (k) is bounded.
From 42 we obtain
& |lm (k+1) - m 20 s x
| (k++1 -b(0)
mm (++1) - m (1) < 0o.
(49)
0=
This implies that m(k) (k) is bounded.
Finally, from 13 we obtain u (k+1) = (k+1) m (k+1) + m (k) and 2 (b) = -b(1) m2))) + m 8-1 Then
1 (k+1) (A) =1b((k++) ((0) m (k) m (k+1) k m (
< [|](k+]) (kk) 1|2 + Il (l++1 - m (x) 113 + | (t) m - 1.
Consoquently, we obtain
| (k+1) - r (1)111 < 20.
(50)
Hence, the sequence (b) (h) m (4) is bounded.
Proof (c): The augmented Lagrangian function (b, u, m) f(b)+g(u) 6 b- - [/ defined
as Ly, R. + (-o0, oo] is proper and lower semi-contimuous, where f(b) 11 xb612 g(u) = ||/|/1+ +|///o
and h(b, u, m) lb mm mm Tf T., u, m) is semi-algebraic, then it satisfies the KT. propreyry
23

---

Published in Transactions on Machine Learning Research (11/2024)
at any point of its domain. Note that both and h are real polynomial functions, which are semi-algebraic functions (Attouch et al. 2013 Bolte et al. 2014p. Both ||ullo and |//|1 have piecewise linear graphs and are therefore remi-algebraic (see Example 3 and 4 in (Bolte et al.
2014). respectively).
Furthermore, consider that 91(2)= = Aa||u//1 and 92 = A(1
lo Their proximal operators have piecewise
linear graphs and are perfectly known objects (Attouch al. 2013 Beck 2017). The proximal operator for gi gi(u)= Aal/uili prox s( (u) [[]] >a]+ sgn(u) (the so-called soft thresholding function) is defined as
u Aa, 0,
2 Aa,
[1c] Aa]+ sgn(s))
if uu Xa,
u + Ao,
u S -
Hence, prox 91() (u) has a piecewise-linear graph and semi-algebraic. The proximal operator for de can be written as
0, if |u| < 22 a),
prox (2(1-a))) = u if |u| > 221(1 a),
if |u| = 221(1 - -a). (0, u),
Clearly, (X(1-a)()) is also piecewise linear and scmi-algebraic. Note that prox (M Hy(u) the so-called hard thresholding operator, is defined as
0,
if |ul < V, if |u| V,
Hy)u) =
{0. u).
if |u| = V,
where v 22X(1-a) a). Here, 91 + 92 is also semi-algebraic.
Consequently, for any nonnegative real numbers À and a, the function s(b) + Aal|ua/1 + A(1 - a.) | lla
Il1b u + mll3 2 [ml] is semi-algebraic. Hence, we conclude that the Lagrangian function in (0 is a KL function.
Since (b (N) m (8) } is bounded, there exists a subsequence {b u m kl converging to stationary point (b",u" m"), where E N. Since the Lagrangian function in 20 KL. function (using the lower semicontinuous property), we have
L,(b", u², m²) < lim (b k .m m)
(11)
1-000
In conclusion, all the conditions (a)-(c) in Theorem 1 hold.
B.3 Proof of Theorem 2
Starting with (PRSM), we update b, m, and a iteratively according to 24a
24d
(k+1) argmin L, (b) (k), ()).,)
(52a)
b
m (k+ = m (k) + b (++1) u (k)
(52b)
u (k+1) argmin L, (k+11) u (2) m (
(52c)
m (k+1) = m (k+ + (k+1) u (k+1)
(52d)
Proof (a): From 52a
since (k+1) minimizes L, (b( (k), m (k) and Lagrangian is strongly convex with
respect to the variable b, 37 holds.
Next, using the augmented Lagrangian function in 20 we compute
Ly (b++1) (k+1) m (.+1) - (b((t))) ( (k+1) m ++ = (b.+t1) - u u(tt))) (m (++1) - m m (53)
24

---

Published in Transactions on Machine Learning Research (11/2024)
Now we rewrite 62d
as
b (k+1) m (k+1) m m+1
(54)
From
52c
and
29
obtain
f(b(+t))) + (b((++1
(b++) + m A b)
(55)
Substituting 52d
into
55
we obtain
41
Again, from 52a) and 29 we obtain
{(((++))) + (b(+++1) w (K) + m (00) 0.
(56)
Substituting 52b
into
56
we obtain
vf(b(k+l)) -
(57)
Using . 54 and
53
bocomes
ml (k+1)) m m =11p) = -11 1 ((b + 1 [( 0.
(58)
Hence,
L (b((et)))) ( ((++))), m (1+1)) - Ly (b(k+.))) (x+1), m + S O.
(59)
Using 20), we have
(''mm" (m*+3 (b(k+)), ukl) +) Ly (1)),, = = (bl++t) - ((()))))((()))) (
ml
(60)
Next, reformulate 52b
(+++1) (x) m++ (8)
(61)
Using 57
and 61
60
becomes
? [| m (11p)) =7 (11 1 ?
Vf(b(k+))) + γ
(62)
Therefore,
(b(++1)), (k) m +++++ L. (b(++1)) (k), m (8)) < [|b(k+1] b 6(4))13)
(63)
where 1 2 0 is a Lipschitz gradient of the function f(b). From 52c we have that
I. (b(+))) u (++1)
m b1 T.y (b (k+1) m k s O.
(64)
Finally, combining
37
59
63 and 644, we get the desired inequality as follows
L, (b(--)) 1) to (k+1) m (k+1)) L, (b(H,, (())),(())) (+)
In (b((+))) (1-1) m (k+)))
T., (b(i+1), (k+1) m b1)
(b(()) u (l+1)
m x - L., (b(+1), u (k), m bb)
(b((+)))) (8) m x 1) L, (b(++1), (k) (4,,m()) (2)
(b((1)))) (k), (m (k)) L., (b()), (x) m (k))
12 [ - b 6001
<
(1111b(++11
1b -b
which is the sufficient decreasing condition
45
25

---

Published in Transactions on Machine Learning Research (11/2024)
Proof (b): The difference lies in some steps to show the boundedness of u and m The rest is the
same as in the proof of Theorem [b). From (54) and 61 we obtain (k+11) b (k+1) m (++1) m R+ and
u (8) ) = b 1))+) m + + m respectively. Then
1| (k++) u (1) = ||m h+ m
(k-1) + m R - m 11
Im (k+1) m
++ 112 || m (10 10)
This inequality can be rewritten using 58 and 62 as
1| u (k)
1 (b(+tt) - 10. -b
Consequently, we obtain
| (k+t) -z (1 < 20.
(65)
Equation 65 implies that (k) is bounded. To show that ml()) is bounded, we analyze the difference m (k+1) m follows m (k+1) m (k) m (k-1) m m+d m+ m Then, obtain
13 S [bb (1b(h+1)-b000 100
ml (k+1) (8)
7
m (M+1) - m (+) xc.
(66)
k=0
Hence, m(() is bounded.
Proof (c): See the proof of Theorem 1 (c).
For the case E (0, 1), all conditions are valid. Hence, for the soquences (b()), (*)), .m ( generated by the SCPRSM scheme 224a 24d), and its Lagrangian given by (20). the three conditions from Theorem [ (a)-(c) hold, achieving a worst-case convergence rate of ( Here, a worst-case O convergence rate indicates that the solution's accuracy, based on specific criteria, improves gradually at rate proportional to one divided by the number of iterations (k) within an iterative algorithm (He al., 2014).
C Datasets
C.1 Simulated QTLMAS 2010 Dataset
The dataset comprises 3226 individuals across 5 generations, including 20 founders (5 males and 15 females) with two observed traits (responses): a quantitative trait and binary trait (Szydlowski & Paczyńska 2011 Each female mates once, producing approximately 30 progeny per birth. SNP data were simulated using a coalescent model on five autosomal chromosomes, each 100 Mbp long. A total of 10031 markers were generated, including 263 monomorphic SNPs and 9768 biallelic SNPs. The continuous quantitative trait is controlled by 9 major QTLs at fixed positions, including two pairs of epistatic genes, maternally imprinted genes, and two additive major genes with phenotypic effects of -3 and 3. The additive genes are positioned at SNP indices 4354 and 5327, whereas the major epistatic locus is at SNP 931. Additionally, dominance locus was positioned at SNP number 9212, with an effect of 5.00 assigned to the heterozygote and 5.01 to the upper homozygote. Moreover, an over-dominance locus was placed at SNP 9404, with an effect of 5.00 assigned to the heterozygote, -0.01 to the lower homozygote, and 0.01 to the upper homozygote. After filtering SNPs with MAF < 0.01, 9723 markers were retained and transformed into onc-hot encoding, resulting 29169 genomic markers. We used the quantitative trait in our study. Generations to (individuals 2326) were used for training, and generation (individuals 2327 to 3226) served test data.
26

---

Published in Transactions on Machine Learning Research (11/2024)
C.2 Real Pig Dataset
The Pig datasel contains data from 3534 individuals, with high-density genolypes and phenolypes for five traits (Cleveland et al. 2012). Using the PorcineSNP60 chip, 52842 SNPss were assessed and filtered to 50282 based on a minor allele frequency threshold of < 0.01. The chosen trait had a heritability 0.58. After adjusting the phenotypic data and excluding individuals with missing data, the final dataset inchuded 3152 individuals and was transformed into one-hot encoding, resulting 150840 genomic markers.
C.3 Real Mice Dataset
This dataset comes from an experiment aimed at identifying and locating quantitative trait loci (QTLs) associated with various complex traits a population of mice. The dataset contains 1814 individuals who were genotyped for 10346 polymorphic markers and two traits: body length (BL) and body mass index (BMT). In this study, used BL trait. After transforming the data into one-hot encoding, the dataset resulted in 31038 genomic markers. This dataset is from the Wellcome Trust and is available in the R. package BGLR (Pérez & de Los Campos 2014).
D Implementation of Baseline Methods
For comparison purposes, we implement the LASSO 3 using the proximal ADMM and SCPRSM schemes (referred to as LASSO-ADMM and LASSO-SCPRSM, respectively) as
((++1) prox (u m (k)),
u (++1) = prox (b (b+1) + (1),,
(67)
m (++1) == (R) + b (k+1)
and
(k) (H), (k++) = prox prox (y (u
m (++ })= (k) +rbb(t+1) (1)),,
(68)
u (k+1) = prox (b+++) + m ( m (k+1) = (k+ + (b(11) u (k+))),
where prox,. 97 (b) = (b)
Similarly, for EN-ADMM and EN-SCPRMM, prox 97 (b) T+P Sy),, where g > a linear combination of the L and 2 penalties (Parikh & Boyd 2013 The closed-form proximal mappings of the SCAD
and MCP 6 penalty functions be found in (Fan & Li
2001 Liao et. al. 2023 Wang & Liu
2024
Yun
let al. 2021
Here, utilize the scaled versions
(b)(b) (a- 11(b) sign(b)aA a
if {b S (1+ if (1 7)} |6| Sa, if [b| > ad,
prox (b) -PrOX (b)
(69)
b
Sxx(b) if b < a7\, m7
ay
prox 9y (b) = prox proxmepp (b)
(70)
otherwise,
with respect to SCAD and MCP, respectively. All iterations of LASSO-ADME LASSO-SCPRSM, EN- ADMM, EN-SCPRSM, SCAD-ADMM, SCAD-SCPRSM, MCP-ADMM and MCP-SCPRSM terminate upon
achieving convergence, defined by the condition (bb**) - u (1) < 3(1 + Il (k) parameter 3 is set to 10-5
|,>,, where the tolerance
27


Reviews:
[{"review_id": "dBTONlXEpK", "reviewer": "Reviewer_uxuY", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\nThe paper is written in a straightforward manner and is easy to follow with minimal prerequisites from optimization. The structure makes intuitive sense and the ideas flow clearly.\n\nIn terms of content, the combination of $\\ell_0$ and $\\ell_1$ penalties is, to the best of my knowledge, novel. The theoretical results are straightforward and the inclusion of different optimization schemes doesn’t come off repetitive. Finally, the usage of the translation functions is interesting.\n\n**Weaknesses**\n\nThere are some weaknesses in the work that can be mainly summarized in the following axes:\n\n- unclear motivation,\n- issues with the theoretical results,\n- limited experiments.\n\n*Unclear motivation*\n\nIt is not clear why we should be interested in using $\\ell_0$ and $\\ell_1$ together. In the introduction some motivation for this is introduced, but it is not convincing. For example, a reason the authors state for using sparsity is to reduce overfitting (first paragraph). In fact, in the second paragraph it is stated that $\\ell_1$ reduces the magnitude of coefficients, and thus reduces overfitting. However, if the desire is to reduce overfitting, a (provably) better choice is $\\ell_2$ regularization. Why is then $\\ell_1$ chosen over $\\ell_2$? Later in the introduction (top of page 2) it is stated that small coefficient values that are produced by LASSO are difficult to interpret and lack information. These statements seem to be conflicting, and moreover it’s not clear why the combination of $\\ell_1$ and $\\ell_0$ would avoid the stated problem of small coefficients. Finally, the elastic net is mentioned (but not compared against). It’s not discussed why the elastic net is not a good solution to this problem, which is rather strange, as it also has a very simple and efficient proximal operator.\n\n*Issues with the theoretical results*\n\nFirst, a comment on the $\\ell_0$ component: in the second paragraph of page 2 a discussion is made about the NP-hardness of $\\ell_0$. Then, at the top of page 3 it is stated that $\\ell_0$ debiases the LASSO and because it follows LASSO then the computation is feasible. Unpacking that further brings up a few questions:\n\n- What does, formally, debiasing mean in this context?\n- Mentioning the NP-hardness of $\\ell_0$ is peculiar to me. First of all, that NP-hardness is never addressed, in the sense of the authors showing how they overcome that hurdle. Second, mentioning NP-hardness seems irrelevant, since such results refer to global optimality. The authors instead use IHT which, if KL or RIP are satisfied, leads to stationary points, not global optima.\n- Finally, it’s not the LASSO that makes $\\ell_0$ “feasible”, it is IHT, and the feasibility itself could be questioned since the proposed methods are 4-20x slower.\n\nIn terms of Theorems 1 and 2, these seem standard proof schemes for Langrangians simply adjusted for the specific Langrangian of (13). There is a technical mistake in Theorems 1d and 2d: the statement says that global convergence is proven, but the actual proofs show convergence to a stationary point. Emphasizing more on this, it is unclear if (12)/(13) are general formulations for any $f,g$ or specifically for (4). The appear to be for (4), but in (19) the derivatives of the Langrangian are taken; at the very least the terminology is highly problematic as the $\\ell_0$ does not admit a derivative, not even a subgradient.\n\nMoving on to the implementation, the introduction of the translation functions is confusing. They are only discussed in A.4, but discussion reads that for affine transformations the proximal is a translation, and thus we can write a translation as $T(\\boldsymbol{u}) = f(\\boldsymbol{u} + \\boldsymbol{m}) - \\boldsymbol{m}$. It’s not obvious to me why this is true, why it is useful, or how it is used in (25). A proof that (24) is equal to (25) seems necessary here.\n\nFinally, the mention of proximal operators in the context of $\\ell_0$ seems misplaced. The authors frequently cite Parikh and Boyd (2013), but that and other monographs (Bach, Jenatton, Mairal, Obozinski, 2011) define proximal operators strictly for convex functions, which $\\ell_0$ is not.\n\n*Marginal experiments*\n\nThe experiments are hard to interpret without significant background in biology. Currently the description of the datasets is extremely technical and it is very challenging to even extract wha the tasks are. Regardless, the performance seems lackluster. In most cases the performance improvement is marginal and comes at a very high cost of 4-20x of increased computation.\n\nMoreover, for LASSO, SCAD, and MCP ADDM is always used. However, in the authors’ methods it seems that SCPRSM is the method that performs more favorably. Given the small differences in performance, it is not convincing that the other methods wouldn’t perform better under similar optimization schemes.\n\nFinally, the main motivation of the paper was to introduce $\\ell_0$ in order to reduce the number of variables in the model, reduce overfitting, and improve variable selection. None of these aims are supported by the current experimental section, which focuses solely on MSE improvements.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a novel regularization scheme that combines $\\ell_0$ and $\\ell_1$ regularization with the goal of retrieving relevant features (through the $\\ell_0$ regularization), while at the same time preventing overfitting (through the shrinkage of $\\ell_1$ regulization). The corresponding optimization problem is defined and then solved through three separate techniques (ADMM, PRSM, and SCPRSM), alongside convergence results. A disucssion about selecting the learning rate and tuning the hyperparameters is presented, and then the models are evaluated against baselines (consisting of LASSO, SCAD, and MCP using ADMM) on a synthetic and two real world genomics datasets.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "There are some changes that would greatly improve clarity. For example ADMM, PRSM, SCPRSM and SCAD and MCP are only introduced in the abstract. Similarly, SNP is only defined in the first few sentences of the introduction. These abbreviations are not used until much later, so when the discussion focuses around these terms it would help a lot in terms of clarity to restate them in full so the readers don’t have to go searching for them.\n\nMore on ADMM, PRSM, SCPRSM, why are you considering these optimization schemes? Howe are they different, what are the relative benefits of each, and why do you use all three? It seems there is no motivation to use them, and even a motivation for (12) (meaning, why was the problem infeasible before, how (12) makes it feasible) is missing.\n\nI mentioned this above as well, but for someone with no significant biology background parsing and understanding even what the task is in the experimental section feels impossible. The descriptions are overly technical which doesn’t seem helpful for the main text and would be more appropriate in an appendix. It would help a lot readability if a simple description or explanation of the datasets and tasks was given and the technical details were moved to an appendix. As a further example of this, in the first experiment an MSE of $\\sim 65$ is reported. In any experimental setting I’m familiar with, such an MSE would indicate a failure on the task. Making the experiments more understandable (and advocating why such a seemingly high MSE is actually acceptable) would go a long way.\n\nYour conclusion states “The WL1L0 method outperforming all known regularization methods”. This is absurd as a statement, as you only test against three of them, not every known regularization method (as I mentioned before, strikingly comparisons against the elastic net are missing). I think it would strengthen your experiments a lot if you included comparisons against the elastic net and unregularized methods. Currently, the emphasis of the experiments is performance. If that is the main purpose, you should test against methods that could potentially outperform yours, otherwise the experiments should be designed in a way that is consistent with your motivation (sparsity, variable selection, overfitting).\n\n**Comments**\n\nFinally, I have some minor comments:\n\n- Why is (29) defined? It’s identical to (27).\n- Have you tried using accelerated versions of optimization, such as FISTA or Nesterov acceleration? These could potentially improve the running times of your methods.\n- This might be a bit subjective, but 3.6 feels overly long to explain a hyperparameter choice.\n- Appendices A, B, C are not named and are all named “Appendix”.\n- What are the assumptions for the Figures 1, 2, and 3 in the appendix? It appears only the norms are plotted, not the error term, so that’s not $F$.\n\n**Mistakes**\n\n- The sentence “In cases where the variables are highly correlated, LASSO may randomly select variables” appears twice in the introduction. The sentence that follows the duplicate sentence is also extremely similar.\n- Zhang (2010) above (9) does not have a softlink.\n- Second paragraph in page 4, two closing quotations ($\\textrm{'oracle property'}$ instead of $\\textrm{`oracle property'}$).\n- Page 9, in the first sentence of 3.6, “relazation”.\n- Above the unnumbered equation before (57): “using the proximal operator for $g_2$ the prox of $g_2$ can be written as”", "broader_impact_concerns": "No broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "zWJB2YhYjj", "reviewer": "Reviewer_rSMH", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. This paper is well-written and easy to follow.\n\n2. The new proposed model can take the benifit of both L1 and L0 norm and have good performance. \n\nWeakness:\n\n1. The problem (4) is derived from problem (2), and I am particularly interested in understanding the relationship between their solutions. Clarifying this connection would enhance the understanding of the model's development.\n\n2. In the experiments, the parameter $\\alpha$ is always small, implying that the L0 regularizer is more influential in ensuring sparsity. It would be valuable to compare the performance of your model against one that uses only the L0 regularizer.\n\n3. Tables 1 and 2 show that ADMM solves problem (14) effectively and efficiently, with a shorter runtime compared to PRSM and a comparable minimum MSE. The inclusion of the parameter $r$ significantly improves performance when $r<1$. I recommend providing a more detailed explanation of the intuition and benefits behind SCPRSM.\n\n4. I suggest introducing the K\\L{} condition before presenting Theorem 1, as some readers might not be familiar with it. This would improve the clarity of the theoretical exposition.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents a new model that combines both L1- and L0-norms (WL1L0) to enhance regularization in high-dimensional data. Despite the model's nonconvex and nonsmooth nature, the authors prove global convergence for several optimization methods, including ADMM, PRSM, and SCPRSM. The model's efficiency and the effectiveness of the proposed algorithms are demonstrated through experimental results.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See weakness.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "cLqwL3aldv", "reviewer": "Reviewer_Vs8U", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: the problem of L1 and L0 learning is important, and the datasets are nicely chosen\n\nWeakness: There are quite a few points that I fail to understand, and the authors should have done a better job explaining them in the paper\n1. I do not understand how the L0 part is being minimized. L0 is NP hard, and how does transforming equation (4) to Equation (12-13) helps minimizing the L0 penalty? This part is not explained. Also, since there is a theoretical complement, perhaps the theorem statements should state in a clear way how the optimization algorithm minimizes the L0 penalty\n\n2. I am not sure about the comparison metric. Is the MSE computed on the test set? Or on the training set?\n\n3. Why do the authors use ADMM to minimize LASSO? There are many optimization algorithms, why ADMM in particular?\n\n4. I think think the authors should also compare with the method in https://arxiv.org/abs/2210.01212, where exact nonproximal algorithm for minimizing LASSO is proposed and may achieve much better performance here\n\n5. The proposed algorithm seems so slow -- taking often more than a magnitude of time to run than any baseline. This appears to be a huge drawback to me", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work proposes a joint weighted L1 and L0 norm regularization using proximal operators and translation mapping technique", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please answer the questions above and include a new comparison", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes a new regularization approach that combines $l_0$ and $l_1$ regularization for preventing overfitting. The paper develops the theoretical convergence guarantee and provides the experimental results with comparison to various baselines. The authors have mostly addressed the reviewers' concerns about the mismatch between theoretical claims and the mathematical proofs, and have fixed errors in certain proof steps. Hence, I recommend to accept the paper.

However, the experimental results indicate that the proposed algorithms incur higher computational costs compared to the baseline methods. Thus, I request the authors to add a discussion about the computational limitations of the proposed method.
Decision: Accept (poster)