Sparse Modal Regression with Mode-Invariant Skew Noise
Authors: Kazuki Koyama
Takayuki Kawashima
Hironori Fujisawa
Venue: Accepted by TMLR ()

Abstract: Sparse regression methods have been widely used in many fields for their statistical effectiveness and high interpretability. However, there are few sparse regression methods with skew noise, although statistical modeling using skewness is becoming more important, e.g., in the medical field. The Azzalini's skew-normal distribution and its extensions are well-used for skew noise. Such skew regression methods have a severe problem with statistical interpretability because they model neither mean, median, nor mode. To overcome this problem, we propose a novel sparse regression method based on mode-invariant skew-normal noise. The regression model is easy to interpret in the proposed method because it always models a mode regardless of skewness. The proposed method is simple to implement and optimize, suggesting it is highly scalable to other machine-learning methods. We also provide theoretical guarantees of the proposed method for the average excess risk and the estimation error. Numerical experiments on artificial and real-world data demonstrate that the proposed method performs significantly better and is more stable than other existing methods for various skew-noise data.

Published in Transactions on Machine Learning Research (09/2024)
Sparse Modal Regression with Mode-Invariant Skew Noise
Kazuki Koyama
kkoyama Dism.ac.jp
The Graduate University for Advanced Studies (SOKENDAI)
Takayuki Kawashima
kawashimaac..iiee ac.jp
Tokyo Institute of Technology
Hironori Fujisawa
fujisawa Rism. ac.jp
The Institute of Statistical Mathematics
The Groduate University for Advanced Studies (SOKENDAI)
RIKEN Center for Aduanced Intelligence Project
Reviewed on OpenReview: https: // openreview.w... net/ /orum! id= 63r6M1JkXm
Abstract
Sparse regression methods have been widely used many fields for their statistical effective- ness and high interpretability. However, there are few sparse regression methods with skew noise, although statistical modeling using skewness is becoming more important, e.g., in the medical field. The Azzalini's skew-normal distribution and its extensions well-used for skew noise. Such skew regression methods have severe problem with statistical inter- pretability because they model neither mean, median, nor mode. To overcome this problem, we propose a novel sparse regression method based on mode-invariant skew-normal noisc. The regression model easy to interpret in the proposed method because it always models a mode regardless of skewness. The proposed method simple to implement and optimize, suggesting it highly scalable other machine-learning methods. We also provide theo- retical guarantees of the proposed method for the average excess risk and the estimation error. Numerical experiments on artificial and real-world data demonstrate that the pro- posed method performs significantly better and is more stable than other existing methods for various skew-noise data.
1 Introduction
Starting with Lasso (ibshiram, (HID)) an l regularization regression method is becoming popular for predictive modeling with relevant features. It has gained significant attention in many fields because of its ability to handle high-dimensional data efficiently and its exhaustive studies of methodological extensions and theoretical properties. (For details, e.g., see Bühlmann & Van De Geen (20T)); (Hasii et all (2M16). As we have more data available and need more reasonable models, an & regularization regression method will likely continue to be more essential for high-dimensional data analysis.
The Lasso is based on the squared error; in other words, it focuses on a symmetric noise. In real-world data analysis, however, we need treat skewness. For example, statistical modeling using skewness is indispensable in the medical field. Hossain & Bevene (20)) explored application of the norw-l distribution in the analysis of microRNA (miRNA) data, which often does not follow a normal distribution. Their findings, derived from both simulations and real miRNA datasot analyses, demonstrated that the skew- normal distribution could enhance the detection of differentially expressed miRRAs.. This enhancement. is particularly noticeable when the data is significantly skewed. Shatter et all (2022) introduced an automated stain normalization framework for histopathology images that uses a mixture of multivariate skew-normal distributions to capture both symmetric and non-symmetric observations. The method can model multi- modal and multiple correlated distributions, regardless of their non-symmetry, and has shown consistent and

---

Published in Transactions on Machine Learning Research (09/2024)
Azzalii (1985)
0.6
0.4
0.0
Fujsown and Abe (2015)
0.8
0-0
2-05 1:5
0:45
0.6
0.2
Figure 1: The difference between the Azzalini's skew-normal distribution (upper row) and the mode-invariant skew-normal distribution (lower row) with the location parameter μ and the different skewness parameter The scale parameter is set to = 1. The blue and black dashed lines in the upper row represent the positions of the mode and in the former distribution, respectively. The mode of the former distribution depends complicatedly on the skewness parameter as well as the location and scale parameters. As a increases, the mode goes first to the right and then back to the left (to the original location). This behavior makes it difficult to interpret the former distribution. The red dashed line in the lower row represents the positions of both the mode and μ in the latter distribution. The mode of the latter distribution is invariant for any As - oo, both asymptotically approach half-normal distributions. The former is stretched vertically, but the latter is stretched horizontally.
superior performance compared to existing methods. How to treat skewness is also topic of interest in other fields, including financial, meteorological, and industrial data (Adrock et all 205 Hao et all 2mg Simola of all, 2MT))
Skew distributions have been studied extensively (Azzalini & Capitanio, 1999
Branco & Dey. 2001 Saln
at 2003 Liseo & Loperhid, 200:: Conzzlez-Pariaa et all 2003 Arellano- Valle & Genton, 2mss Arelanno
Valle Azzalini, 2006). However, despite the intense need for skewness in real-world data analysis, there currently exist few relevant & regularization regression methods. The most significant study is Chen al (201E), who assumed regression model for the location parameter of the well-used Azzalini's skew-normal distribution (Azzalimi, (185] and then proposed an & regularization regression method with an efficient parameter estimation algorithm. Even if we move away from the & regularization, there exist few regu- larization methods; e.g., Can all (2122) recently proposed the ridge regularization regression method for the mode the Azzalini's skew-normal distribution. This lack of attention may be because we would gen- erally deal with skewed data using non-linear pre-processing with a power transform such as the Box-Cox transformation (Box &r. Cox,, (ire)) and Yeo-Johnson transformation &
2000). However, these non-lincar transformations may not be suitable for regression model with skewed noise because such trans- formations provide symmetric distributions of outcomes but not symmetric noises due to asymmetric feature distributions. Furthermore, non-linear transformations can make it difficult to interpret regression model since the non-linearity changes the meaning of an observation unit.
In this paper, we focus on a regression model whose noise is assumed to follow the mode-invarian, skew- normal distribution introduced by Fujisawa & Abd (201), and then we propose the regression estimator of the mode (location parameter) that is defined the minimizer of the negative log-likelihood plus the L regularization, aiming to achieve both skewness modeling and statistical interpretability. The mode- invariant skew-normal distribution is subclass of the Aaansfrrmation-of-Saal (ToS) distribution (Jones) 2014 2016) and is consistent with a normal distribution when not skewed. As the name implics, its mode is invariant regardless of skewness, where its location parameter coincides with its mode. Therefore, the proposed method always models a mode regardless of skewness, making easy to interpret the statistical role of sparse features. This is oquivalent to the reasonable assumption that a noiscless situation has the highest probability. Despite its helpful properties, there are few studies with the mode-invariant skew-normal distribution besides the basic maximum likelihood estimation (Fujisawa & Abe 201)
2

---

Published in Transactions on Machine Learning Research (09/2024)
The proposed method differs significantly from Chen et (201z))'ss 's regression model for the location param- eter of the Azzalini's skew-normal distribution. Figure M illustrates different shapes between the Awzalini's skew-normal distribution and the mode-invariant skew-normal distribution. Note that the location parame- ter of the Azzaliní's skew-normal distribution is neither mean, median, nor mode, which mcans that the role of Chen all (201)) regression model is uncertain and difficult to interpret.
Our main contributions are as follows:
Highlighting the attractiveness of the mode-invariant skew-normal distribution and assuming noise to follow it, we propose a novel regularization regression method having both skewness modeling and statistical interpretability.
Revealing partial convexity, we indicate that its optimization and implementation are relatively straightforward, although the proposed method involves complicated non-linear transformations.
The proposed method theoretically guaranteed terms of its excess risk and estimation error by non-asymptotic analysis under mild assumptions usually adopted in the theoretical analysis of Lasso-type problems.
This paper is organized follows. In Section the proposed method is introduced. Tn Section the theoretical guarantees of the proposed method are provided by non-asymptotic analysis in terms of the excess risk and the estimation error. Such theoretical studies for skew regression models have not been shown so far. In Section M. numerical experiments on artificial and real-world data are demonstrated. The results obtained that study indicate that the proposed method can serve as a viable alternative to Lasso. In Section 5 concluding remarks are given.
2 Proposed Method
2.1 Skew Distribution
We use the mode-invariant skew-normal distribution (Fujisawa & Abc 2015 with the parameter vector Φ = o,, E R+ IR:
1
)
fly
xb) ((ul()=--
(1)
where @ is the probability density function of the standard normal distribution, Ta special transformation function of scale, described later, and μ,, o, and the paramcters of location, scale, and skewness, respectively. The function To must satisfy
(a( > 0, La (v) + (-v)= =
(2)
where yo the inverse function of ra. Note that the condition (2) ensures a monolone-inereasing property of (u). The model (M) is unimodal with the mode-invariant property Mode [f(y] +)] μ for any and
Let Ha(v) == qa (v) - While Ha (or ra) can be freely designed if it satisfies the above conditions, this paper employs the following special function Ha because of its better analytical treatment:
2
(o(v) :
(4)
a
where Pa : exp(-a?).
Remark 1. The above function Pa slightly differs from that proposed in Funsawa Es Abe (2018) because the Fisher information matria non-singular at a 0 when we use the above function but singular when we use the function proposed in Fujisawa Abe (emr)). The non-singularity is usually necessary obtain some theorelical properlies. A related proposilion oppears in Proposition
3

---

Published in Transactions on Machine Learning Research (09/2024)
2.2 Problem Description
We consider the linear regression model with the skew-unimodal noise, given by
y = X g +e,
(5)
where y E R an output, X E is a P-dimensional feature vector, E R.P is a parameter vector, and e follows a skew-unimodal distribution f(e 00,0,a) a). Let (xX) X B,o,a]" The distribution of y can expressed by fly 4b (X)),, and the mode X
Let the log-likelihood function be denoted by 1(0 X y) log ( (x)) with the parameter vector 0 [8 "0',,,a) a.aa] Let DN == (Wwn Xn) be the iidd. data. We define the loss function by
N Xn)
((0) DN) =
Σ
1(0
N
n=1
1
Σ
log [Vn (x)))
(6)
N
Since we are interested in sparse regression, this paper focuses on the Lasso-type problem
min A(0) (2) 111
(7)
where À 2 0 is the regularization paramcter.
2.3 Optimization
(i) is non-convex for but convex for when and are fixed, as described in the following theorem.
The loss function
Theorem 1. The loss function E(e Dx) is convez with Lipschitz continuous gradient for B. In particulor, if SN = N - X, is positive definite, ((O) (D)) strongly conver for B.
The proof is given in Appendix R. Using this theorem, we develop a problem-specifie optimization method based the block coordinate descent method. This method allows separate the optimization problem into convex and non-convex and take advantage of the convenient property described in Theorem M.
The update algorithm is summarized in Algorithm M We update 8, a, and alternately. First, we fix (a,a) and optimize B. We employ the inntorization-MMnmmmii with acceleration techniques (Jamsiddaa & Jemnnan 1997 Sun et all 2016). The derivation of this MM algorithm given in Appendix C As shown line of Algorithm I. we can rewrite the B update as another Lasso-type problem and then well-known software, e.g., the sklearn. linear model package of Python. Although
(MM) algorithm using the Taylor expansion
Next, we fix (B, a) and optimize o, and finally fix (8, and optimize the non-convexity for and a, their optimization problems are at most one variable. Hence, use trustable software, c.g. the scipy, optimize package of Python. In this paper, we employ the L-BFGS algorithm (Liu & Nocedal, (188)), which an iterative method for solving non-linear optimization problems. In particular, for the inequality constraint > 0, we can utilize the L-BFGS-B algorithm (Byrd et all [oss Zhu et all (997)) which extends the L-BFGS algorithm handle bounded constraints.
need to consider easily
Let {(0) {(0) Dx) + \||B||1- Let dom ? be the effective domain of i. Following the definition of Tseng (200T)), say that 0 E dom is stationary point of & if for any d E RP+2
i(0 + ad) {(0) 0.
lim inf aj0
Algorithm M guarantees convergence to stationary point, as described the following theorem.
Theorem 2. Let (0(1)
sequence of 0 after updating in Algorithm D. Then, every cluster point (0(")))_,1..... I,.,, is a stationary point of F.
4

---

Published in Transactions on Machine Learning Research (09/2024)
Algorithm 1 Optimization of the proposed method
Require: prper-parameter À and initialized B3, o,
1: while until convergence do
2: while until convergence (of (3) do
3: for n ,,...,N .N do
4: Zn - (yn X 3)
+ Pa
Wn 22 + Pa
< X B+ a
in
1tpa
Un +
W/r
7:
end for
8: B - argmin (1-(1-12) Σns (in x, 3) + A|1
1+
8 end while
9:
10: - argmin E(0 (D)) with L-BFGS-B (or other valid algorithm)
a:
11: + argmin !(0 Dx) with L-BFGS (or other valid algorithm)
12: end while
The proof is given in Appendix D. Here, we describe how to set initial values in this paper. The regression parameter B started from the estimated coefficients of the Lasso model for the same data. Using the residuals of its model, we initialized with their sample standard deviation and a with -1 or 1 corresponding to the sign of their sample skewness.
2.4 Related Work
The most famous skew distribution is the Azzalini's skew-normal distribution (Azzalim
implemented in current standard scientific computing software. Many skew distributions were developed using the Azza- lini's basic idea (Azzalini & Capitanid, 1999 Branco & Dey, 2001 Sahu et all 2013 Liseo Loperfido, 20133
Fanzález=Fariass all, 2004 Arellano- Valle Genton, 2008 Arellano- Valle & Azzalini, 2000). The closest method to the proposed method Chen al (2012)) which assumees a regression model for the location parameter of the Azzalini's skew-normal distribution. However, the location parameter in these skew distri- butions is neither mean, median, nor mode; the mean, median, and mode have complicated combinations of all the location, scale, and skewness parameters. Even if the location is modcled by X / for the Azzalini's skew-normal distribution, it is difficult to understand the role of features. The proposed modeling makes it easy to understand the role of features via the modeling of mode. The proposed modeling suitable for an interpretable one in a high-dimensional setting, but the modelings using the Azzalini's skew-normal distribution are not.
3 Theoretical Results
This section provides theoretical guarantees for the proposed method by non-asymptotic analysis. First, we present some basic properties of the first, second, and third derivatives É(*b | y) = log fly | (b) Next, using them, we show that the excess risk has quadratic margin. Finally, we show the upper bounds of the average excess risk and the estimation error, and then we verify the favorable properties of the proposed method. All proofs of propositions and theorems in this section given in Appendices B Gl and H.
3.1 Basic Properties of Likelihood
We present some basic properties of the first, second, and third derivatives of b(1b v). These properties are the basis of non-asymptotic analysis.
5

---

Published in Transactions on Machine Learning Research (09/2024)
Suppose that X' [X] is included in a bounded space of RRP. We assume that the parameter spaces v and Q (0) are bounded in the following sense. Hereafter, we suppose that the following assumption is satisfied.
{2)}
Assumption 1. For some positive constant K, we suppose
V c v {bb l S K, 1/K S |o| S K, |a| {K).
(9)
θ c θ =0 (x)) E L for any X E X}.
(10)
The true parometers sp and 0 are included in v and θ respectively. Let the score function (first derivative) of 1(ab v) log fly ab) be denoted by
a
(D = Oub
l(ab lv).
(11)
Proposition 1. We have
| S C12lul C1 + Ci.0.
(12)
where C2 C1,1, Cr,0
positive constonts.
The score functions of regression parameters μ and are bounded by the first and second-order polynomials of ly)- This is the same property as the normal linear regression model. We may expect that the score function of skew parameter a is bounded by a third-order polynomial of lyl because the skewness usually depends on the third moment, However, it is bounded by a second-order polynomial of ly]. Proposition II implies that the proposed model easy to treat a as well as a.
Let the Fisher information matrix of fu( (1) b)) be denoted by
2²
I(tb) = Est
Oub l(tb (l(b v)
(13)
Let x [XT go, 00, Ta!T. The Fisher information matrix [(()(()) is uniformly positive definite for any X, in other words, inf Amin (I(bb(X)) 0, where Amin (A) is the smallest eigenvolue of A.
Proposition
If the model is simple, we can easily prove a similar proposition to Proposition 2 because the Fisher in- formation matrix is easily calculated. However, the model (M) is not simple, and it is not easy to verify Proposition 2, as seen in devised proof of Appendix [a We found that the constant factor fa proposed by Fujisawa & Abe (2 does not satisfy Proposition at a = This is why a constant factor Pa: is slightly modified, as described Remark M, and then we can prove Proposition D
Proposition 3. There ezists some function G3(y) such that
g
sup (61.2..))6(1.23333
Oy,, Oyh, Brh I(ab v) s Gs(v).
C3 C3:=sup sup / (a(l))(v) (x)) <0. x
(14)
3.2 Excess Risk and Margin
Let the Kullback-Leible divergence between foo and fvr which is called the excess risk, be denoted by
E(bb 26) = Esyo [og for (1) (y) [v)
(15)
Proposition 4. For any > 0, there caists some constant d. such that
inf ||
inf (X)laae
E(bb (x)) 2 8..
(16)

---

Published in Transactions on Machine Learning Research (09/2024)
Proposition Hl shows the identifiability condition, which is usually necessary to estimate the true parameter by a statistical method correctly. The identifiability condition is often assumed implicitly, but it is precisely verified here because the model (I) is complicated due to skew noise.
Theorem 3. There crists some constant C > O such that
& (ab 25 (X)) X || (x) 2 C
inf
(17)
Theorem & says that the proposed method has a quadratic margin, implying the error & (ib (x)))) presents a sharper behavior around (x) than the well-used squared error || - (5(()) up to a constant factor.
3.3 Convergence Rate
The following is well-used condition in non-asymptotic analysis for the & penalized method, which is called the restricted eigenvalue condition.
With L > 1, define the So) So) -restricct Assumption 2. Let So == {p: 52 # 0} and So == [I,...,P)|.., PM...
eigenvalue constant by
S
8² (L) = 18ggllish 18 lI |881
inf
(18)
We set == 6 and assume 0 < So) 1.
The squared loss of the simple linear regression model for low-dimensional case is strongly convex be- cause the Hessian matrix of the squared loss concerning the regression paramcter vector 3 is Sy and it is positive definite. This property provides the strength of the least squared estimator. However, for the high- dimensional case N < P, Sy may not be of full rank; in other words, its smallest eigenvalue may be zero, and then the squared loss not strongly convex. The non-asymptotic analysis reveals that if the smallest eigenvalue of Sy under the restricted condition || || L||8s, 112 is positive, such as Assumption 2, we can
show theoretical guarantees of the & regularized estimator (Bihhlman & Van De Geen 2014 Hastie et all,
2016
In Assumption 2, we set I. to simplify later notations. We replace it with any constant, greater than one by adjusting how to choose the regularization parameter. Based on this assumption, we have the following important theorem.
Theorem 4. Suppose that Assumption is satisfied. Fir X1, XN. Let the average excess risk and the empirical process of the loss function D be denoted by
E(O 00)
Σ E (v()))) | (xx))) n=1
(19)
N
VN (0) = (xb),,) vnn) - E (+(((),)))) b)))
(20)
Let |0]|, : 1|11]1 + ||o, a]T 1l.. Fir y 2 1 and Ao > 0. Let
J VN (0) - Vn (09)
S 5
(21)
| 00 v X
If the regularization parameler X is chosen to satisfy X 2~λo, then we have on T
E(θ + 2 1811, 7
(22)
where C is the constant defined in Theorem - = Ri, So), and so #So-

---

Published in Transactions on Machine Learning Research (09/2024)
Suppose
Ao 0 (nog N2
log N P)
(23)
This order is necessary to prove Theorem 5 and is relevant to the following two corollaries. From Theorem EI. we can obtain convergence rates regarding the average excess risk E( (*) and the estimation error of B.
Corollary Suppose À = 2 with 22 and the same assumplions as in Theorem + Then, we have
$oollo 10g W
E(θ 0d) OF
(24)
N
(og((V PJ)
|| 1 OF
soooog N)?
(25)
N
The orders of the average excess risk [241) and the estimation error [25] are (log N) and (log N)² larger than those of the ordinary Lasso, respoctively. These are the same as some well-known complicated models, such as the linear mixed effects models with regularization (Schelldorler el. all, 20T1). A similar result is also seen in Corollary
3.4 Feature Selection
When the non-zero coefficients are sufficiently large, which the so-called beta-min condition (Biihaam 2012,, we can show that the estimated set of non-zero coefficients == {p 3p + 0} includes the true set of non-zero coefficients So with high probability.
Corollary Suppose 2y with (28) and the same assumplions as Theorem g. Suppose the non-zero
coefficients of Bo are sufficiently large so that
min 10 > O ( so(log PeSo
log(N V PJ)
(26)
N
Them, with high probabilily, il holds $ So-
As in Corollary M. the order in 20 is (log NY² larger when compared to the ordinary Lasso.
3.5 Reasonable Success Probability
The statements after Section S9S hold when the event, J in 241 occurs. Finally, we show that the event. J occurs with high probability.
Theorem 5. Suppose À = 20 with (2x) and the same assumptions as in Theorem 4 The event J occurs
with probabilily al least - 6, where with some positive conslands C4,1, C4,2 > 0,
(1
log log(N V V)
1 N
8 C4,1 exp
+
(27)
C2
4
Experimental Results
This section demonstrates the performance of the proposed method with synthetic and real-world datasets.
4.1 Skew-Noises with Mode Zero
The simulation model xT + e. The feature vector X,, E IR P=110 was generated randomly with
each element following the continuous uniform distribution 1, IR. The true parameter was set

---

Published in Transactions on Machine Learning Research (09/2024)
2.0
X10
1.0
Sparity
190
Nodel
3.5
2-3
0.6
12-
Propoced
Shant
Leac
Lamse=1)
No
M65 05
MS6 VK
Sporsiy 211
100
Model Sice (ti
F3 Scone
Propoaed
Leaci11)
3-
0.6
0.0
x19"
x16-1
0.-
0.5-
3.5
12-
0.6-
0.8
10-
Preposes
Char
No.44
VDR IV
x10
VS2 I
Spars 21
100
Model Sux (N
F1 Score
1.5
Chen
8
Lamce
Trepose Chena Leoe U03017)
Propeued Chens Laoz Loso
Prosssed Duni Laoz Leso-1
Propzued Chans Leso Leocs-1)
08
10
Figure 2: Comparison of four methods with four different types of simulation noises in terms of five evaluation measures. From top to bottom, each row results from Data 1, Data 2, Data 3, and Data 4, respectively. The first column plots the probability density function of the noise. The second through fifth columns are box plots of some scores. The last column is the mean F1 score for each threshold (horizontal axis). In all figures, the proposed method, Chen+ (Chen all 2014), Lasso, and Lasso with the Yeo-Johnson transformation are rod, green, bluc, and yellow, respoctively. All experiments were conducted for 50 runs with different random seeds.
80 = [1.0, -0.9, -.,1, ... T) E 100 with 10% non-zero coefficients of varying magnitudes and signs. Let the random variable following (M) be denoted by SN o, a). Let the log-normal variable be generated from four types of distribution: SN 1, (Data 1). denoted by In,) o). The noise exp(-0.53) (Data 3), and In (n, 1) exp(-12) (Data 4). Data 3 SO 1,2) (Data 2), (n.V(0,
and Data 4 are the log-normal variable adjusted to mode zero. These are misspecified models. The noise distributions are depicted in the first column of Figure 2
We compared the proposed method with the ordinary Lasso (Tibshiram (996)) (*Lasso"), the ordinary Lasso applied the Yoo-Johnson transformed output ("Lasso+YJ") and another skew-noise Lasso (Chen ct all, 201E) ("Chen+") that assumes regression model for the location parameter of the Azzalini's skew-normal distribution. The tuning power parameter in the Yeo-Johnson transformation was determined by maximum likelihood ostimation. The sample size was set to N 500. We conducted 50 experiments with different random seeds. For each trial, the regularization coefficient À was adjusted by 5-fold cross-validation based on the log-likelihhoe loss, in which the numbers of training and validation ion data for each trial were set. to 400 and 100, respectively.
Five evaluation measures are shown in the second through six columns of Figure D. The second column is the Mean Squared Error (MSE) for the prediction: MSF(O) Entest (On 2 where the MSE

---

Published in Transactions on Machine Learning Research (09/2024)
MSE (yl: N
MSE (yl: P
MSE (B): N
MSE (8): N
Proposed Lasso
Proposed Lasso
10"1
Proposed Lasso
Proposed
3.8
Lasso
10-2
3.6
3.4
10-2
2017
3.0
100
100
200
300
100
500
600
800
s0n
1000
10
100
100
200
300
400
5000
600
100
800
900
1000
Figure 3: Box plots of MSEs for Data 2 with P = 100 fixed and various sample sizes N. The vertical axis is a logarithmic scale. The horizontal axis is the sample size N. All experiments were conducted for 100 runs with different random seeds.
estimated from 2,000 test data generated under the same conditions. The third column is the MSE of the estimated coefficients: MSE($) == |13- Since "Lasso+ YJ" is based on the Yeo-Johnson transformed output, the resulting estimator 3 does not generally converge to 8 Therefore, MSE(3) is expected to be inferior to other methods. The fourth column is a measure for sparsity by Hurley & Rickard (2004) as the most robust. measure founded on a weighted sum of estimated all coefficients for evaluating sparsity. The measure range is [0, 1] E IR; as a coefficient vector is sparser, its value is higher (see Appendix D for details). The fifth column the size of the estimated model, which is defined the number of non-zero coefficients:
Model (ize( := # 13, # 0 1.15p SSSP P}. The last column is the F1 score (Van Rijsbergen, 2018)) defined by
Precision Recall
F1 Score := 2 Precision + Recall
(28)
where Precision and Recall are defined by TP/(TP FP) and TP/(TP FN), respectively, and TP, FP, and FN mean the true positive, false positive, and false negative, respectively. The positivity (negativity) of the estimated coefficient is defined by 1 > ôr for a threshold ör [0, 1] (horizontal axis). As the feature selection is more successful at each threshold, the F1 score is larger. This means that better method presents upper curve.
As a result of Figure D. the proposed method outperformed the comparative methods in all the cases, indicating better prediction and feature selection with a smaller number of features. We see that Chen+ is inferior to Lasso in Data 3. This might be because the Azzalini's skew-normal distribution could not adequately model such skewed noisc.
Figure 5 shows box plots of MSEs for 100 trials regarding Data 2 with P 100 fixed and various sample sizes N. The regularization coefficient À was tuned by 5-fold vols-valid with 80% training and 20% validation data. The Lasso is a benchmark against the proposed method. The proposed method presented stable behaviors for N < P as well as N > P.
In Appendix 12, we also show the results with a linear regression model without regularization, assuming the Azzalini's skrwanormal distribution, skew-l distribution, and Cauchy distribution for noise (Arellaume
V Azzalim 2013 Azzalim &r Arellano- Valle, 2013 Azzalim & Salem 2020), on Data 1 through Data
4.2 Normal and Azzalini's Skew-Noises
We explore both normal and Azzalini's skew-normal noises under the same conditions as in Section In Let the random variable following the Azzalini's skew-normal be denoted by SN SN (j4, a.). The noise was generated from additional four types of distribution: N/(0, 1) (Data 5), (0, (Data 6), S0, 1,4) (Data 7), and (0, 1,6) (Data 8). Notably, the Azzllini's skew-normal noise is not
10

---

Published in Transactions on Machine Learning Research (09/2024)
2.0
1.30
x10
1.0
Sparity
190
Nodel
3.5
2-3
1.4
0.6
12-
Shant
1.06
Leac:
Lamse=1)
Nr.c
MSA
VSC IV:
Sporsit
100
Model Sic Wii
F3 Score
5.0
Propoaed
Leac-1)
2-
0.6
0.8
5.0
+19
x16-1
2.0
0.5-
2.5
52-
0.6-
0.0
10-
0.4
Preposes
Char
No.e
102 19
x10
V52 I7
Spars 21
100
Model Sux (1
Score
Chan
1.0
Lamce-1]
Trepose Chera Leoe UB9917)
Propeoed Dhens Laoz Leso
Progzeed Dun Laoz Leso-1
Propzued Chan Leso Leocs
08
1.0
Figure 4: Comparison of four methods with four different types of simulation noises in terms of five evaluation measures. From top to bottom, each row results from Data 5, Data 6, Data 7, and Data 8, respectively. The first column plots the probability density function of the noise. The second through fifth columns are box plots of some scores. The last column is the mean F1 score for each threshold (horizontal axis). In all figures, the proposed method, Chen+ (Chen al 2014,, Lasso, and Lasso with the Yeo-Johnson transformation are rod, green, bluc, and ycllow, respoctively. All experiments were conducted for 50 runs with different random seeds.
mode-invariant ie., its location parameter does not coincide with its mode, which significantly violates the assumptions of the proposed method.
The results are shown Figure In Data 6, Data 7, and Data 8, the proposed method outperformed the other methods with the difference becoming more significant Chen+ was slightly worse. This would be because the objective function of Chen+ is hard to optimize due to troublesome non-convexity with complicated combinations of all parameters. Similarly, for Data 3 (skew data), Chen+ was worse than Lasso. In comparison, the proposed method can present stable behaviors because the objective function convex for many variables B E R.,. although is non-convex for only two variables a, a R, as shown Theorem 1.
Data 5, the four methods presented similar behaviors, as expected.
increased. Surprisingly, Chen+ presented similar results to Lasso, although
4.3 Real-World Medical Data
We applied the proposed method the following two medical datasets: PDGFR (Platelet Derived Growth Factor Receptor) consists of N 79 samples and P 320 features, where the outcome the ability
inhibit. PDGFR phosphorylation (Guha Jure 2ME MTP (MelTing Point) includes N 274 samples
11

---

Published in Transactions on Machine Learning Research (09/2024)
Table Results of PDGFR (N 79, P 320).
MSE
Sparsity(33)
Size
Proposed
7.00 10-1 (3.16 10"1) "
9.90 10-1 (3.16 8)
5.30 10 (1.84 10%)
Chen+
10.9 x 10-1 (2.54 x 10
8.18 X 10 (3.72 10-1)
8.17 x 10 (1.82 x 10")
Lasso
7.28 x 10-1 (4.79 x 10
9.85 x 10 00.2211000
x
8.90 x 10 (5.57 x 10")
Lasso + YJ
12.7 x 10-1 (4.22
9.85 x 10-1 (9.82 1
8.93 x 10 (5.57 x 10")
Table 2: Results of MTP (N = 274, P 1142).
MS)
Sparsity(3))
Size (8)
Proposed
8.43 10-1 (1.65 1)
9.91 10 (2.05 ")
2.09
10a
(5.45 10")
Chen+
10.8 x 10-1 (1.88 x 10"))
9.72 x 10-1 (4.58 x 10-2)
5.93 x 10 (7.47 x 10")
Lasso
23.4 x 10-1 (81.9 10
9.82 x 10 (6.00 10
4.25 x 10
(13.9 10")
Lasso + YJ
11.6 x 10-1
9.82 x 10-1 (6.66 x 10
4.25 x 10
(13.9 x 10")
and 1142 features, where the outcome is the melting point of drug-like compounds (Karthikeyan el all, 2015). These datasets have many features compared their sample size and may hold skewed noise, described later. We compared the proposed method with "Chen+" (Chen ct all, 2014), "Lasso" (Tibshirani, (199)), and "Lasso+ YJ", in the previous section.
We used 80% of the samples for training and the remaining 20% for testing. Then, we tuned each regular- ization coefficient with 5-fold cross-validation using 20% of the training samples (i.c., 16% of all samples) as validation data. These samples were generated randomly, and 30 trials were conducted with different random scods. We calculated MSE(Ú), Sparsity of B, and Model Size( and evaluated their performance.
Table M and Table Z show the results of PDGER and MTP, respectively, with the means and standard devi- ations (in parentheses) of the evaluation measures. For both datasets, the proposed method simultaneously achieved the smallest prediction error and the smallest model size. In particular, for MTP, the proposed method reduced the model size to less than half when compared to the other methods. Regarding the resid- uals of the test data after training with Lasso, the mean of the normalized skewness of the residuals was -0.81 for PDGFR and 0.75 for MTP.
In Appendix K we also applied the proposed method to more easily interpretable data, specifically the Engineering Graduate Salary (EGS) prediction data. (Aggarwal et all 20m).. EGS has fewer features than the two datasets above, and the meanings of all features are specifically provided. These properties allow us to consider the validity of the estimated active features. The result shows that the proposed method outperformed the comparative methods and could select more reasonable features.
5 Conclusion
In this paper, we have proposed a novel sparse regression model whose noise is assumed to follow the mode- invariant skew-normal distribution. The proposed method always models a mode regardless of skewness and highly adaptable and interpretable to skew noise. We have also shown that the proposed method is straightforward to implement and optimize, and it has thcorctical guarantecs for the average excess risk and the estimation error by non-asymplotic analysis. The results of numerical experiments demonstrated that the proposed method achieved both skewness modeling and statistical interpretability simultaneously, even for high-dimensional data, highlighting its significant effectiveness.
As a next challenge, we focus on exlending to the mode-inwariant skew-l distribution. To our knowledge, no study has explicitly dealt with this distribution besides its basic idea being mentioned in Fujsawa & Abd (201). The first problem would be to prove the positive definiteness of the Fisher information matrix, whose proof could be complicated, such as in the proof of Proposition R Its optimization algorithm could also be problematic. However, these may be overcome in terms of t distribution being scale mixture of
12

---

Published in Transactions on Machine Learning Research (09/2024)
normal distribution. It is also interesting to extend from univariate to multivariate noise for some machine learning algorithms. Fortunately, multivariate mode-invariant, skew-normal distribution has already been proposed by Abe & Fujisawa (2011)) suggesting that we could quickly tackle and formulate models.
Acknowledgments
The authors would like to thank anonymous reviewers for their insightful comments and suggestions. This work was partially supported by JSPS KAKENHI Grant Numbers JP17K00065 and JP22K17859.
References
Toshihiro Abe and Hironori Fujisawa. Multivariate skew distributions with invariane through the transformation of scalc. Japanese Journal of Statistics and Data Science, 2:529-544, 2019.
Christopher Adcock, Martin Eling, and Nicola Loperfido. Skewed distributions in finance and actuarial science: review. The European Journal of Finance, 21(18114):2533-255 2015.
Varun Aggarwal, Shashank Srikant, and Harsh Nisar. Ameo 2015: A dataset comprising amcat test scores, biodata. details and employment outcomes of job seekers. In Proceedings of the 3rd IKDD Conference on Data Science, 2016, pp. 1-2, 2016.
Reinaldo B Arellano- Valle and Adelchi Azzalini. On the unification families of skew-normal distributions. Scandinavian Journal of Statistics, 33(3):561155) -6454-444 2006.
Reinaldo B Arellano- Vallan and Adelchi Azzalini. The centred parameterization and related quantities of the skew-t distribution. Journal of Multivariate Analysis, 113:73 90, 2013.
Reinaldo B Arellano- Valle and Mare G Genton. On fundamental skew distributions. Journal of Multivariate Analysis, 96(1):93 116, 2005.
Adelchi Azzalini. A class of distributions which includes the normal ones. Scandinavian journal of statistics, pp. 171 178, 1985.
Adelchi Azzalini and Reinaldo Arellano- Valle. Maximum penalized likelihood estimation for skew-normal and skew-t distributions. Journal of slatistionl planning and inference, 143(2):419 433, 2013.
Adelchi Azzalini and Antonella Capitanio. Statistical applications of the multivariate skew normal distribu- tion. Journal of the Royal Statislical Sociely: Series B (Stalistical Methodology), 61 (3):5779 602, 1999.
Adelchi Azzalini and Mahdi Salehi. Some computational aspects of maximum likelihood estimation of the skew-t distribution. Compulational and Methodological Stalistics and Biostalisties: Contemporary Essays in Advancement, pp. 3 28, 2020.
George E.P Box and David R Cox. An analysis of transformations. Journal of the Royal Statisticnl Society Series B: Statistical Methodology, 26(2):211 243, 1964.
Márcía D Branco and Dipak K Dey. A general class of multivariate skew-elliptiaa distributions. Journal of Multivoriate Analysis, 79(1):99 113, 2001.
Peter Bühhlann. Statistical significance in high-dimensional linear models. In 8th World Congress in Probability and Statistics. 8th World Congress in Probability and Statistics, 2012.
Peter Büihlaann and Sara Van De Geer. Statistics for high-dimenssioal data: methods, theory and applica- lions. Springer Science & Business Media, 2011.
Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu. A limited memory algorithm for bound constrained optimization. STAM Journal scientific compuling, 16(5):1190-122)) 1995.
Xingyun Cao, Danlu Wang, and Liucang Wu. Performance of ridge estimator in skew-nooe normal mode regression Blaiisies-Simmlaaiion and Compulation, 52(3):11644 1177, 2023.
model. Commenications in Statislics- Simedation
13

---

Published in Transactions on Machine Learning Research (09/2024)
Lianfu Chen, Mohsen Pourahmadi, and Mehdi Maadooliat. Regularized multivariate regression models with skew-t error distributions. Journal of Statistical Planning and Inference, 149:125 139, 2014.
Hironori Fujisawa and Toshihiro Abe. A family of skew distributions with mode-invariance through trans- formation of scale. Statistical Methodology, 25:89 98, 2015.
Graciela González-Parias, Armando Maminnuez-MMoinne and Arjun K Gupta. The closed skew-normal distribution. Skew-elliptical distributions and their applications: journey beyond normality, pp. 25 42, 2004.
Rajarshi Guha and Peter Jurs. Development of linear, ensemble, and nonlinear models for the prediction and interpretation of the biological activity of set of pdgfr inhibitors. Journal of chemical in information and computer sciences, 44(6):2179-2D)) 2189, 2004.
Songhua Hao, Jun Yang, and Christophe Berenguer. Degradation analysis based on an extended inverse gaussian process model with skew-normal random effects and measurement errors. Reliability Engineering & System Safety, 189:261-270, 2019.
Trevor Hastie, Robert Tibshirani, and Martin Wainwright. Stalistionl learning with sporsily: the lasso and generalizations. CRC press, 2015.
Ahmed Hossain and Joseph Beyene. Application of skew-normal distribution for detecting differential ex- pression to microrna data. Journal of Applied Statistics, 42(3):477 491, 2015.
Niall Hurley and Scott Rickard. Comparing measures of sparsity. IEEE Transactions on Information Theory, 55(10):4723-4741, 2009.
Mortaza Jamshidian and Robert I Jennrich. Acceleration of the em algorithm by using quasi-newton methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 59(3):569-55) 587- 1997.
MC Jones.
Generating distributions by transformation of scale. Statistica Sinica, pp. 749-771, 2014.
MC Jones. On bivariate transformation of scale distributions. Communications in Statistics- Statistics-Theory and Methods, 45() 588, 2016.
Muthukumarasamy Karthikeyan, Robert Glen, and Andreas Bender. General melting point prediction based on a diverse compound data set and artificial neural networks. Journal of chemical information and modeling, 15881998,, 2005.
Bruncro Lisco and Nicola Loperfido. A bayesian interpretation of the multivariate skew-normal distribution. Statislies & probabilily letters, (4): 401, 2003.
Dong C Liu and Jorge Nocedal. On the limited memory blgs method for large scale optimization. Mothe- matical programming, 45(1):503-55) 528, 1989.
Julien Mairal. Optimization with first-order surrogate functions. In International Conference on Machine Learning, pp. 783-797 PMLR, 2013.
Sujit K Sahu, Dipak K Dey, and Márcia D Branco. A new class of multivariate skew distributions with applications to bayesian regression models. Canadian Journal of Statistics, (1(2):229-150, 2003.
Jürg Schelldorfer, Peter Bühlmann, and SARA VAN DE GEER. Estimation for high-dinensional linear mixod-effects models using l--penalization. Scandinavian Journal of Statistics, 38(2):197-214, 2011.
Sobhan Shafeei, Amir Safarpoor, Ahad Jamalizadeh, and Hamid R Tizhoosh. Class-agnostic weighted nor- malization of staining in histopathology images using spatially constrained mixture model. IEEE trans- actions on medical imaging, (9(11):3355 2020.
Umberto Simola, Xavier Dumusque, and Jessi Cisewski-Kehe. Mcasuring precise radial velocities and cross- correlation function line-profile variations using skew normal density. Astronomy Astrophysics, 622: A A131, 2019.
14

---

Published in Transactions on Machine Learning Research (09/2024)
Ying Sun, Prabhu Babu, and Daniel P Palomar. Majorization-miniiizznn algorithms in signal processing, communications, and machine learning. IEEF Transactions on Signal Processing, 65(3):794 816, 2016.
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267-288, 1996.
Paul Tseng. Convergence of block coordinate descent method for nondifferentiable minimization. Journal of optimization theory and applications, 109:475-494, 2001.
Cornelis Joost Van Rijsbergen. The geometry of information retrieval. Cambridge University Press, 2004.
In-Kwon Yeo and Richard A Johnson. A new family of power transformations to improve normality symmetry. Biometrika, 87):9)954 959, 2000.
Ciyou Zhu, Richard H Byrd, Peihuang Lu, and Jorge Nocedal. Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization. ACM Transactions mathematical software (TOMS), 23 (4):550 (40:550-660 1997.
A Upper Bounds on (a(v)) and |ra (u)
Proposition 5. We have (Ha(v) < palu), [va(v)] < Pa. [.(e)l < Pala), and [e)] < 3pa From < 1, we also have [H)( < [u], [h)[v)] 1, [v)( < lad,, and (x((v) < 3a².
Proof. It holds after simple calculation that
d Ha Hv( = Pa
ha(v) de
(29)
v1 a²n?²
Pal
(v( du
Ha(v)
(30)
(1 + 22
d
a
Ha (v) (1 + 222)
(31)
du3
We clearly see [ha(v)] < Pa and (,(v)) < Pala). Since Ha(v)= = Jod (u)dd,, we have (H) < Palo). We also have
1 (1 a1+a202) (1+ a2222 S 3p
lowl
[M) 3p
(32)
Proposition 6. Suppose that Assumption M holds. Then, r'd (w), (u), and (") are bounded and [ra(u)] are bounded by clul with some constant c.
Proof. For a = 0, since (u) = u, this proposition clearly holds. Hereafter, we consider the case a * 0. By differentiating u = (r(u)) with respect to u, we have = d'a (u)) and then
1
(a(w)) (r(ra(u)) 1 (ra(u))
(33)
From Proposition F. we have
1
0 I+pa Pa
<r((w) S 1 Pa:
(34)
15

---

Published in Transactions on Machine Learning Research (09/2024)
We can easily see that 1/2 Pa with some constant 1 because Pa 1 (1/2) exp(-a²) and is bounded from Assumption M This implies that (u) is bounded. From (0) = 0, we have ra Jo (s)dss and then
Ira(u)) - 1 Pa
|ll < clul.
(35)
By differentiating (8% with respect to u, we have
N/ (ra(()))) (u)
(u)
(36)
{1+
By differentiating the above with respect to u, have
(u)= (un===-"(r((u))rra... 3h, (uraa
(37)
From Proposition and the boundedness of a, can casily see (u) and (" are bounded.
0
& Bez Proposition 7. Fir a + 0. For the definition in (+)) the upper bounds of Ha()))
and dos Ha (t) can be described by linear function of |u]. They can also be described by linear function of da Ha(w)),
8a
8² 8a (v) and 8² de Ha( can be described by quadratic function of [o], constant, and quartic function of [w], respectively. They also be described by quadratic function of lu), a constant, and a quartic function of lu), respectively.
ul/. For 0, the upper bounds of @
Proof. Analytically, for fixed a + O,
Ho (v) = [aexp(-aa)
1
+ avI 1 + az
Ha(v)
da
Pa
< (eexp(=a2)
+ Ha(u),
(38)
Pa
8²
(pa (pa- 1 2a²)
exp(-ca? -(ae)ee))(((a)
2
da Ho (v)
Ha(v)
pa
(1+a2
a² + a
eexp(--or
Ha(v)
Pa
av1 + a2²²
Pa 2a eeppp=aa
2v3
Ha(v)
9a a²
aeepp_aa'
1 1
Ha Ha(v),
(39)
Pa
a
16

---

Published in Transactions on Machine Learning Research (09/2024)
03
(- 20(ps 22)))ppp exp(-a²)
20(Po exp( exp(-2a")
20x(1 + +a)) exx((pxX) Ha(w)
Ha(
0a3
Pa
Pā
Paa
30v
8²
+
Ha(v)
(1 + 22))
a1l + 22)2)
a / + 221
(Pa 2a²) expppxa*
(
(acxp((aa)
1
3
Pa
(1 (1+a2
a2v1+
Pa
av + a²
Ha(v)
eexp(-or aa
Ha(v)
Pa
av1 a²u²
20 expl (-a²)
+ 20
20 expl- -202) +
20(1 + Pa) exp(- (v(
Pa
Paa
Pi
1
48
2v3
+
+2
Ha(o)
a
125
9
Pa 2a exp(p -2
2v3
[eeep(-220) 222))
1 Ha(v)
+3
+
9a +
+ a
Pa
a=
Pa
(aexp(-a*
a²)
1
(v).).
(40)
Po
a
Therefore, since 82 Daz Ho (v) and z Ho (v) can be described by a linear function of lol. Furthermore, from 55 their upper bounds can also be described by linear function of lu].
[Ha(v)] < Pa |u| from Proposition it holds that the upper bounds of 88 H. (v)
Next, consider the case O. Let g(x) vI x² 1. Letting the m-th order derivative of g(x) be denoted by (x²), we have
g (x) V1 + x
(41)
(2) (z)
(42)
(1 + x2)
3r
)x) (z)
(43)
(1 + x2)
Then, it holds after simple calculations that
glav)
lim
(44)
① g(av)
(_ g(vv)
+ v (aw)
lim
lim
a da
a-10
a²
1
2
(45)
8 g(av)
[2(av)
2u ()
v² (cv)
lim
lim a0
+
a
da²
a²
0
(46)
8
g(av)
6g(av)
Gv g (av)
302 g (2) (av) (av)
83 g (av)
lim
lim a a0
a-+0
3 a
3
a
2
a
3
(47)
4
17

---

Published in Transactions on Machine Learning Research (09/2024)
It also holds that
lim
Pa
lim
1
1
expl
a-
a-0
2
1
(48)
2
0
lim
Pa =
lim a (xepp(-a²)
a-30
da
a-0
(49)
82
lim
Pa
lim (1 exp(xx))
a a-0 da?
a a-0
=1
(50)
lim
Da
Po
lim a-1
2a(2a² - 3) exp(-xa
(51)
Note that (z( Pa
g(ar)
Using from
1 through
5
we have
g(av)
g(av)
lim
Ha (v) = lim
Po
+ Pa
a-40 da
a
da
a
da
1
2
(52)
8
8²
g(av)
g(av)
0² glav)
a-0
lim
da
Ha(
lim a-0
da dazPa-
Pa
+
da
Pa
da a
+ Pa
da² a
(53)
glou)
02
á g(au)
8² g(av)
g3 g(av)
lim a-0
3o
(v)
= lim
a-0
da
Pa
+ Da Pa
da
+3 da Pa-
da
+ Pa
da
3
2
(54)
2
Hv(v) as a - 0 can be described by a quadratic function of |o|, a constant, and a quartic function of |ol, respectively. From (8F), their upper bounds can also be described by a quadratic function of lu], a constant, and a quartic function of [u]. respectively.
Therefore, the upper bounds of à Ha (a))),
8² Ha(w)), and da
0a
daz
Proposition 8. Fix a + 0. For the definition in (B)), the upper bounds of da (u)), de ra and das rasaau) can be described by a linear function of |u], quadratic function of ju], and cubic function of a
0² ul, respectively. For a = 0, the upper bounds of da (u)),, dae ra)). and des ra(u)) can be described by a quadratic function of [u], a quartic function of [u], and a seatic function of |u], respectively.
Proof. Using the Eulcr's chain rule,
i (u) = - θ
Ha(
(55)
du
da
8²
Ø2 - (u) Ha(v), Da?
0a? Ta(u) = -')(
Halv))
(56)
Da
03 0n (
( a da
② Ha(v)
8
= -- (u)
Ha(()) (r) Da
Ha(
r'((u)
H. (v).
(57)
Oa.
Therefore, from Proposition E and Proposition the proof is complete.
18

---

Published in Transactions on Machine Learning Research (09/2024)
B Proof of Theorem 1
Let (u a) = (u)). Let La(a u) := log f- a). For 0, since the optimization problem can be regarded as an ordinary least squares method, this theorem clearly holds. Therefore, in the following, we consider the case + 0. From (y 1 and μ XT B, since
0 0303T e(0) Dx) =- 1 N
Σ
11,0, 0,
Xas Vn)
0308T
n=1
N
02
x
(a un).
(58)
No²
Ou?
it is sufficient for Theorem @ verify the convexity of la(2 un)) for un Letting en := (un),, have
0²
1 02
(a un)
2 dug r (un)
Oug
(vn)) + To (un)rr(un)))))
+ ha (Un) Un h1, (vm)
(59)
(1 ho (vn)))
where the last equality holds from the calculation in the proof of Proposition Since (vn)l < Pa from Proposition E, the denominator of 5 is bounded by
0 < (1 - pa)² < (1 (a())) < (1 + Pa) < 8.
(60)
The numerator of is obviously 1 for Un Consider the case en # O. The numerator of is expressed
ha (Un)
+ ha(v)) vnlla 1 Un
Un
Cn
1 - Pa
(61)
a
and then
0 < Pa
<1 (a())) (,,,((n)
< Pa 2.
(62)
Combining 54
and
h2
have
1 1- Pa
1 1 ha(un) - Unn)(un)
1 Po [1 - Pa)3
0
< 0G.
(63)
+ Pa)
(1 ha (Vn))
Thus, there exist some constants Cr, == (1--aa/(1+pall 0 and Co == (1 pat/p)/(1- pa) 3 > 0, satisfying
N CL
2²
Cy No² n=1
T (()) DN) Σ X X
Σ X
(64)
2 n=1
From Assumption 0, a is bounded and then we see that Cr. is bounded above zero and Cy is bounded below. 2²
In particular, if Sy is positive definite, the minimum cigenvalue of agast {(O Dxx) is positive from (i Hence, {(0) (Dw) is strongly convex for B. The proof is complete.
19

---

Published in Transactions on Machine Learning Research (09/2024)
C Derivation of MM Algorithm
Consider fixing (a, a) and optimizing B. According to Fujiswaa & Abc (2015), the function ra corresponding to (
is
(au + Pa Pa v (an Pa) + +(1 p%)
(a + 0) (ax = 0)
(u = a(1-13)
(65)
11
Let (5(8) = ((e) (N) A|ll|1a with fixed (a, a). For a + O, the optimize problem (+) rewritten as
arg min (F( arg min
log fyy 31 A||
N
2
arg min
Zn
Pa Wn
11131111
2NV
a<11
p&)
==1
N
22 "
20 wn + (2 + - P%a) 1||
arg min
2N
a² P2)2
mm1
N
22) - 2p (1+
arg min
*118111
2N
a²(1 Paj2
12
1 Pa
2pa
arg min
2n Wn
A|1
(66)
2Na2
p22
1
P
where (3) = _ (Yn - XT B) + Pa and (B) = V +
Pa. Thus, there exists some constant Cp.
satisfying
1 P2
N
2pa
(x(B)= 2Naa pPa)
2n
2n wn
+ +11BB111 CF.
(67)
1
p
==1
Note that 2n 2n
Xx and
0g
wn
En Xn, Letting (z(3) = 2/ o wr
1+ 2n Wn, we have
09
0
(9)
2a
2n X, +
0(1 22)
2apa
Wn
Wp
Z
2
X,
2a
2apa
Pa
2n X +
210,
X,
(68)
0(1 P2)
Wn
0303T T (g)
2
2a 2
X
XT
2a Pa:
2a²pa
2
1 Pa
2n
X X
(69)
σ²
02 + 12)
w2
Wn
Since
1-p2 (2
2n um
12 as
Zn - 100 and
8 Izn
2
1 Pa
Zn
3(1 22)2
> 0,
(70)
w2
Wn
w
1-p2 (2+
we have
En w.
2. Then, using 0 <
1/2 from 0 OSpa Pa 1, we have from
I+p0
g2
402
0
in (3) x
x X,
(71)
0303
σ²
20

---

Published in Transactions on Machine Learning Research (09/2024)
From (in and (Z)), we construct surrogate function La(B ((4) of (x(3) near g based the
second-order Taylor expansion as follows (see Lemma 12 of Sun el all (2016) for details):
1 po (r(B) - 2Na p2)2
Σ
N
(3((0))+ +
Sin (8-310)
+ o² (8 X,,,, -8-1033000 - 33(0)
20
x,
03 B=g(t)
n=1
+ 1118111 + CF
= 16 (()).
(72)
Therefore, from and T
we can update g00) at the f-th itcration as follows:
arg min (a(B ( g(++1) arg
= arg min
+ Pd
N
oi,
(g - 3(0) +
20
o2 (B - ) ( (B
x,,X
g(0)) -
+ >11l11>
B
2Na 1
P
08
n=1
B=g(t)
arg min
1 Pa
N
Inn x 2
+ 1||B11
(73)
No
-Pay
3
n=1
where On = X (+)) +
2a
2n
pr 1+p2
Wn
+
2012)
18y
with (1) = 2n (Bit)) and (1) == w (B(1)).
For a == 0, the optimize problem (d) is rewritten as the ordinary Lasso, i.e.,
N
arg min
Σ
(yn- x] + 1113111.
(74)
2No2
==
Proof of Theorem
First, we fix (ar a)) and consider the convergence the MM algorithm for B. We will prove that the surrogate
function (a (8()) in Appendix D is first-order surrogate of (x(3) near ((0)) based on Definition 2.1 of
Mairal 20%
From
22
Lo( (()) is obviously majorant function of ey ( From h.
and
we have
82 0333T
Cy No²
N T
Cv 0²
(r( x
X,X
Ammxx I,
(75)
z:LF
8² 0338T (c( (8())) 2(1 Pa)
2(1 + P:.) 02 - 22)2 Amax
X
(76)
No² - 22 2
ziLo
where Amax(Sx) (55)) is the largest cigenvalue of SN : 1/2 X Hence, 8 (F and 83 La ((t)) are
Lipiscchttzz continuous and LG- Lipschitt continuous, respectively. Let ( g(t)) == la( ((0)) - (g).
Then, o (B g()) is also (LF + La)- Liisshiizz continuous because for any 8x and By:
θ L (Bx ((()) 09
0g (n |00)
<
08 ea [8(0)= _
① 03 ( 180)
cr (r((3xx
l (8x)
03
< Lx Byll + L Byll
(Lp + Lc) Byll.
(77)
Moreover, we have en (p(0)) (g0)) 0 and
83 n (B(1)
2
(g()) 0.
21

---

Published in Transactions on Machine Learning Research (09/2024)
Therefore, (B (()) is first-order surrogate of (x(3) near ( From Proposition 2.2 of (Marral (2013) we have for any 2 1 and some positive constant R
o, + LG)
+
(78)
t+2
where {(0) := E(0 Dx)++ 1 and 3* is a minimizer of & with fixed (a, a).. Hence, it holds that g - 8 00.
Next, we regard Algorithm D the block coordinate descent (BCD) method for the three blocks (8,a, a) and consider its convergence based on Theorem 4.1 (c) of Tseng (2001). When using the cyclic rule (see Section 2 of Tseng (2пл1]) to update the BCD method for B blocks, Theorem 4.1 (c) of [seng (2000T states the following.
Theorem 4.1 (c) of Tseng (2001). Assume that the level set O (0 i(O Dx) S ((0 (x)) is compact and that & is continuous on @o. 1f i/E........ On) has at most one minimum in Os for b= 2..., ,..,., B 1, is a coordinatewise minimuum point of f. In addition, if is regular at 0t, then Ot is a stationary point of L.
then every cluster point 0T of [0
We assume that there always exists (Mn, X,) E DN such that Yn x 3 > O after updating 8 by the MM algorithm, and also, there always exists (ym, x E Dx such that yn. x, g < 0. Under this assumption, it holds that 8(0) Dx) - do as a a a oc, a + -0, or a +
Let dom& l and dom be the effective domain of and É, respectively. We in dom & and the level set @ is compact. From Lemma 3.1 of [Teng (20mT)) 8 is regular at any O dom since dom l is open and is Câteauuxdiiferrettaable on dom f. Therefore, from Theorem 4.1 (c) of Tseng
that {(0 Dx) is continuous
(20mI)) if ē(e Dx) has at most one minimum in a, then Theorem Z holds. We will prove that this necessary condition holds.
Let x(u) == ur. urr(())a(u)) We have
((u) ra (w)ra (r)) (r)((()) = rr (w)r(()) + u (r's (u)) + (ra(u)r(())
(79)
Since (u) > 0 from
and (0) = 0 from
39
it holds that
sgm(ra(u)r (w)) = sgn(r(u)) = =gg)
(80)
Since (ra(4))) ra(u)r"(u) > from iy and it holds that
sgn (u ((r(()))== -rou)))) ra(w))((w)) sgn(u).
(81)
Hence, we have
sgn(x'(u)) = sgn(u).
(82)
Let Un := 9r X For yn - _ > O, Un E (0, o0) decreases monotonically as a increases. Since x(uu,) > 0 for un, (0, 00) from (62),, x(u)) strictly decrcasing function for a. For yn X, < 0, un (-20,0)
increases monotonically as a increascs. Since x' x((u,) < 0 for un (-oc, 0) from (S2) x(u,) is also a strictly
decreasing function for For Un X
x(un) = (0) 0.
22

---

Published in Transactions on Machine Learning Research (09/2024)
Thus, do e DN) is strictly increasing function for o, because
0 DN)
1
i
((0
1(0
X yn)
do
N
da
un (Un (un))
No
n=1
Σx E((mn)
N
N
(83)
No
n=l
Morcover, there cxists only at= 5) satisfying E x(un) N, because
lim
Σ ((un)=
lim Σ x(un) + lim Σ x(un)
010
Un
Un
20,
(84)
N
lim
x(un) = lim
Σ
x(un)+ lim
Σ x(un)
00000
10
Mn
n=1
On
Un 30
++0
(85)
á E(O DN) = and the only minimum of ( | Dx). The proof is ā is also the only solution satisfying do complete.
Note that from Theorem M. if Sy is positive definite, f(O Dv) is strongly convex for and has only one minimum in B. In this case, be optimized first by swapping lines to 9 and line 10 in Algorithm and then every cluster point in sequence of 0 after updating is stationary point.
E Proof of Propositions
E.1 Preliminary Propositions
Proposition 9. The M-th order moment of U ~ fa(u at) is finile. The M-th order moment of Y ~ f(y X TB,,,a) a) is also finite.
Proof. From Proposition 5. it holds that
By-, [U = n" (ra(t))
/ (u
ha (v))do
$ / |
$ (10|+ [4
< (l0l + (1+
(1 + pa) m / lo)
= (1 + pa p -
M oc.
(86)
23

---

Published in Transactions on Machine Learning Research (09/2024)
Furthermore, since y ou + μr. it obviously holds that
Byws [Y]
EU-s. [[oU x ()
m o* (x (3 Th
RUns. [U]] < co.
(87)
k
k=
Proposition 10. Suppose that Assumption 1 is satisfied. Then, the third-order derivative of the oog-likeldhood 8²*
function 1(0) X, y), more precisely, 30,, 00.2 00:3 by a polynomial of lyl-
w) 1(0 x,y)
for any (i,,i,,))) E .I,...,PP P+ are bounded
Proof. The scores of the log- likelihood function 1(0) X, y)
0 x, du du 0g
0 du 0μ L. u)
1(0 "0
03
ra(u)r' ( (u)X,
(88)
0 du 1(0 x,9)
(a [u) + Oo
log
0o
du do
ra(u)r' (u)
(89)
0
8
1(0 x,9) 9)
(a u)
da
da
0
(u) (u).
(90)
da
The second-order derivatives of and with respect 0 consist of
2² da To (u), daz Ta(), das (a()), oa (r)(
a dr))),
0
u, (u), ra (u), (u)) ("(())
da Daz (u).
(91)
as terms involving u. From Proposition B and Proposition E, the second to eighth functions of (M are bounded by some sextic function of u. Consider the remaining last three functions of (II). We have
2 da r (w) (r" Hv) da θ
(92)
0
da ("))) (mm)(u) Ha(v),
(93)
da
0²
i Hal())
2
82
0a ra (u) (u)
(1) Ha(
(94)
du
da2
From Proposition E and Proposition 8 the absolutes of the above three functions bounded by quintic function of lu]. Since u (y u)/a and Assumption M, the proof is complete.
Proposition 11. Let E(ga(v) = (7)). Then, have
[u- Ev~o(ob(()(( (V))
(95)
Proof. On the basis of integration by substitution,
Eu~So(la)][[(U) = ((w))(ra(u)ddu
f (
Ev + h. (V)))
(96)
7
24

---

M Published in Transactions on Machine Learning Research (09/2024)
For condition A, let 1[A) be the indicator function, which means 1(A) if A is true and 1{A) 0 if A is false.
Proposition 12. For ~ (N0, 1) and a sufficiently large positive constant M, it holds that
exp
./2 (Mm-1 + O (M
(-1²)
(m 0)
B Th 1 { > MH -
(97)
exp (- T
(m = 1, 2, 3,
Proof. For m 0,
E ( M)] =2 / ((v)do
JM
2 f
Φ(v
M)du
Ju
M
5 2 exxp
ó(v)do
M2
exp
(98)
For m 1,2,3, 4,
E [ M}]
vo(v)dd
M
V
2
M
S exp
(99)
2
E
5[v 1 > M}
=2 2
8² (')(v)du
M
2
M
M exp
[
([V]
>
MH
2 M exp
M
M
exp
2
2
M2
- (M + 1) exp
(100)
1
MF]
=2
"6(()du
2
M²
M exp
+ 2E [[[]
> MH
π
2 M exp
M
M
<
+ 2 exp
2
2
M
$ (M2 2) exp
(101)
2
25

---

Published in Transactions on Machine Learning Research (09/2024)
E B[IV'' 1 [[V]
Mi] =2
/ Vo(())
M
V Ma
2
exp
M
+3E [ {1V>> > M]]
2
2 M exp
M2
(- A/2))
$ V
+ +3(M + 1) exp
(-2 2
2
M2²
S (M + 3M + 3) exp
(102)
2
Proposition 13. For U ~ fa(u] a) and a sufficiently large positive constant M, it holds that
2exp
(- Ma 5
(m = 0)
B Th 1 { MH S
(103)
4 (M + O 22) exp
(-
M² 8
(m 1,2,3,4) 4)
Proof. From Proposition E. we have [v]] S |o| + [Ha(w) < 2/μl and |1 + ha(w)) < 1 + [ha(v)] < 2. In
addition to this, using Proposition I and Proposition D2. it holds that
E M Th (UU) > My = E (laa((V)) m (1 + ha (V)) flga (V)I > MY]
M 2
< 2mm E [ 1 |V >
ccx
M2
(m 0)
2m+1
m-1 2
M
+0 ((4)
exp (- 8
(m = 3,4,3
2exp
M²
(m = 0)
8
(104)
(M + O (M) exp (-
M² 8
(m = 1, 2,4)
Proposition 14. For Y ~ f(y () and a sufficiently large positive constant M, it holds that
{
2 2exp (-& ( (xx (M= )²)
(m 0)
E [ Th 1 (I > M
(105)
4K (M +0 (2) exp (-& ( )2)
(m- = 2,
where K is based on the parameter space (
Proof. Note that |x| S K and a < K from
Using Proposition IK it holds that
E ([ 1 ( > M2 = E || + oU m
[12 + oU > M)]
S E [K + K (U
[K + K | > MY
|U { |0| >
M K
mm
C.E.E.
K
( (1(M- M- - KK) 2
( K
< 2Kk exp
m + 44Km
(M/ K
)
M K (-2) )
1 11 K) (
m C.
+0
exp (-8
K
K
2 eexp
(- ( (M* 2)
(m 0)
(106)
4K (Mm-1 0 (M (-1 (
(Mx (2)
(m 1,2,
26

---

Published in Transactions on Machine Learning Research (09/2024)
Proposition 15. Let the upper bound of (I22) be denoted by Gs(v)- Let
Gly M) == Gs(y)1 (Go(y) > M}.
(107)
Then, for fly () and a sufficiently large positive constant M, there exist some positive constants Cs.1, Cs and C5,3 such that
[GY| M)] & Cs CxM(C
(108)
E (Y M)] s C3 exp(- exp(-Cg,3M).
(109)
Proof- From Proposition MI
know
Gs(v) C1 + Ci,1vul C1.0 > M,
(110)
and then we have
C1. + VR, + 4 C1,0)
|
2 ss,.
(111)
20,2
where C3,4 is some positive constant. From Proposition
E [GY M)] < E CCuM)
| VCo,uu]]
40 +0()) exp
1
VC3 K
K
S SCs exp(-Cs,MM,
(112)
where C5,1 and C5 some positive constants.
E (V M)] < E [C12ll1] + Ci.lul+ Cr.0 ( {1V1> VC
< K ((VCsal))* +0(a)) exp
VCo.aM --M K
K
s C -G,3M)
(113)
where C5 is some positive constant.
E
Proposition 16. Let G(y) be the function defined by (100] Let
Fly M) = G(y M) + E [G(Y M) X].
(114)
Let My
3 log N for N exp Co,3
27C2.2 6
V
3
ee. Then, it holds that
C5
2 N
F(Y, MN)
(115)
270 C
3C2
Proof. From N 2 exp
V
Ve, we have
Cs.3
3
Cs.2
v3Co.1
< Vlog N < 22
(116)
VCs.aN2
C².
27

---

Published in Transactions on Machine Learning Research (09/2024)
Therefore, using Proposition 15 it holds
loog N
2 2og N
F(Yn MN)
PP
(G(Y, MN) E(G( MN Xx)]))
N
N
N
1 N
2log N
P
G(Y, MN)
Σ
E
E Mx) X,)
N
N
N
2
2log
P
G(Y, MN)
C5,1 M
Cxp(-Cs,,MN)
N
N
Vlo N
V3Cs,1
P
G( MN)
2 2vlog No
N
N
VC
log N N
P
G(Yn
MNN
N
S E
(V MN)
log
N
2
I Co exp( (Cpp())
log .
33C3,2
1
N NologN
(117)
E.2
Proof of Proposition
Proof.
Note that u = (y 3)/o.. This proposition is derived immediately from Assumption II.
SS
SS
( Proposition D. and Proposition R
E.3 Proof of Proposition 2
Proof. For simplicity of notation, let w(v) := de (c).. From the scores
S
and
the Fisher
information matrix can be described as follows by using Proposition II
E
82
1(bb
13
G2 VVn(())
1
V2Ha
V2
V2H2 + V1
V2Ha
V2w
H
1
0
0
0 0
+ (V)
-V2w
V2H
2
0
0
V/2
V2Ha
-V2W
1
Ev-o(e)
V2H,
V2 +3 +3
-V²Haw
20
-V2w
V2H,,w
V2
0
V2
V2Ha
-V2w
1
EVao(())
V2H
V2H2 +1
H22 V2W2
202
-V
-V²Haw
1
Ha
1
E BVroRo(1)
Ha. H2 ++
-Haw w²
(118)
20
- -Haw
Hereafter, since (26(v) is a probability density function on E R, the random variable is regarded
V ro(()) for notational simplicity. Note that matrix (ITS) is clearly positive semidefinite, only need
28

---

Published in Transactions on Machine Learning Research (09/2024)
to check that its determinant is positive. For - 0, we have
a
lim Ha(v) = lim Pa a-0
a-0
v1
22 +1
= 0,
(119)
v²
lim (()) = σ lim exp(-a>) a-+
(V1+a222-1) + + VI+
a-0
(v++a=+++1)
σ
(120)
4
and then we see that the determinant
(188) positive because
Ha H2 +1
w - -How w²
1)-
1
0
i BV [V2
det EVNoRd(e)
Ha
0
1
O
w
- Haw
2 (o [V2]
0 16 Ev-o2((1) [V4]
2
(Evwuig(o)
[V]]
E
~d)(() [ve]
16
3
0²
> 0.
(121)
Consider the + 0. We have
1
Ha
))
E Hal
[
det EvErso(o)
Ha
H2 +1
-Haw w²
=
E [H
E [H2]
+
E. H
-
Haw
W
-E [Hw]]
E W
(E H
E [w2]
E H
E
[w]]
E
H E w
+ 2E [Haw]
[H. E [
E
H
+ (E
[
E (w])
E [Ha] |H [w] 2E
H
H
B|w]
E²
[Haw] + V [w]
= (E -E
- E 1 (()) + [w]
V [w] 2ov2 [Ha.w],
(122)
E² (Ha]) sign of (122), we evaluate V [H.] concretely, starting with (M). i.e.,
where the incquality is based on HHH V [w] = [E [H3]
EE [W]]
E² [w]) > 0. To examine the
1 12
Vyneld(s) [Hta] P2 V vool((v)
V ab
v22 11
1
2
BVmoap(e)
V2+ a
R22 ~oo()())
1
0
3
E
EV
V
V
(123)
a²
29

---

Published in Transactions on Machine Learning Research (09/2024)
[V V2 + +
Here, we focus on the behavior of Ev
with respect to a. We have
0
1 a²
a
+ 2²
da a²
V
1
BvnoRd(e)
da JR
1 + 2012
u² P(()dd
0 JR 80
+ a²
v² (uv)do
1 2nx2
a
20²
b(v)de
(1 + 2a?)
a v1 + a?
+1+ + a²
2
a²
a
(1 + 20
V1+ a2
(1+ a2
22
a
(v)dd
(124)
(1 20²)
(1 + 22)))
where the fourth equality is derived by
| (())du = La vI 22 (o(v)) du
3.
(R)
/
302
a² (1 + 0222)
= 2 lim
(V1+0210
Φlo)do.
(125)
R+000
v1+ R2
Hence, 8a 2+ Ve (oo d(v)
0
1
Ey
V V2 +
is monotonically increasing function of az, Moreover,
a²
1
1 / V22
1 E2
v a²
lim
2 +
2 + ad(e)
1 V
a²
a²
1.
(126)
Thus, we have
E2 v
VV²+
1 a² >
1
(127)
a²
and Vono=o(o) [Ha] < <a < from 12 For a + 0, Hv(v) and w(v) aren't obviously constant func-
tions. Then, V VywoPo(e) (V) < 0o because the upper bound of H. is a linear function of lu] from Proposition @ Therefore, since 0 < [Ha] < and < [(w) < co, it holds that
Cov2 [Ha, < Ha W W WWw
(128)
From 122 and 1223) we see that I($) EE
is positive definite for any sp. Consider the case X = (X 30 o a)). Since (X) belongs to the bounded set V, we have inf Amin (I (x))) >
1(ab Y)
+b
0.
E.4 Proof of Proposition 3
Proof. From Proposition with some natural number M and some positive constant Cam for m =
0.1...., M, we have
g3
M
max (1
1(0
X,y)
Camlly!
(129)
Oyb Oeh, Oei.
m==
Since the m-th order moment of Y f(v X A, a) is finite from Proposition B. the proof is complete. 7
30

---

Published in Transactions on Machine Learning Research (09/2024)
E.5 Proof of Proposition 4
Proof. Applying the Taylor expansion around the point wp for (TE), there exists some parameler veclor ab E {tep" + (1 t)pb E I2 E [0, 1]} such that
ae
0 #2 5/(b +6)TT
-
E(bb 2ba) E(1p") | 22)0) Tes
(eb 260)
(eb 26").
(130)
10 10
Excess Risk 1
satisfies E(60)) 25) = and
arp (260)
08
200) = 0 by its definition. Because I(xb) is positive
definite from Proposition 0 its smallest cigenvalue Amiu (I(z)) is positive. Thus, the quadratic term of (L.C)) is evaluated as
0.
( - bb)
(+b - Φ6) 2 Amin (I( - v
(131)
OrpārpT
Therefore, letting C2 = inf Amin (I(b(X))) 0, based on
inf XEx #E0:110000 1600:110~6611222 1230
inf
E(bb 1/))) 2
inf
inf
Amin
-
Xex $EB.1100
2²
inf Amin (I({)(X))
xex
C²
(132)
2
can take 8, A'O, 0 as &e satisfying Proposition !
F Proof of Theorem 3
If Condition 1, 2, and 3 in Section 9.4.1.1 of Bühhmann & Van De Geen (201T) holds, the proposition holds immediately from Lemma 9.1 of Bühhmann & Van De. (2 These Condition 1, 2, and 3 correspond to Proposition G 2 and Hl which were already shown, respectively. The proof complete.
G Proof of Theorem 4
Tf Condition 1, 2, and 3 Section 9.4.1.1 and Condition in Section 9.4.2 of Bithlman & Van De Geen (201T) holds, the proposition holds immediately from Theorem 9.1 of Bihhlann Van De Geei (2 These Condition 2, and correspond to Proposition G. A and EI, which were already shown, respectively. Condition 4 correspond to Assumption Ø The proof is complete.
H Proof of Theorem 5
Theorem P is similar to Corollary 9.1 of Büihlmann & Van De Geed (2 The difference is model. The former model is regression model with the mode-invariant skew-normal noise, which belongs to a class of generalized linear models, but the latter model is a finite mixture model. The proof of Corollary 9.1 of Bihhlann & Van De Geed (2 is based on two parts; one a model-specific part, and the other is not related to model. The model-spocific part consists of Proposition 9.4 and Lemma 9.3 of Bühlmann & Van De Geed (2ITT). We prove these propositions for generalized linear models. Using them, can apply the (2TT to prove Theorem F.
same proof technique as in Büblmaan Az Van De Geen
Proposition 9.4 and Lemma 9.3 of BBihlmann Xr Van De Geer (201 show how to obtain an upper bound of the score function and how to set Mx. For our sparse regression model with the mode-invariant skew- normal noise, the upper bound of the score function obtained in Proposition M. which corresponds to Proposition 9.4 of Bühhmann
Van De Goed (201I). In this case, following the approach described Bihmaan & Van De Geen (2ITTT) we set My Ollo N)) and then we can derive Proposition MR which corresponds to Lemma 9.3 of Bithmann & Van De Geen (20TT). The proof is complete.
31

---

Published in Transactions on Machine Learning Research (09/2024)
Additional Experimental Details and Results
1.1 Definition of Sparsity
In Section El we used the measure for sparsity (denoted by "Sparsity"), as introduced by Hurley & Rickard (2009). In this section, we briefly review its definition.
Given a vector B E R²", we sort the coefficients from the smallest to largest, i.e., B(1) S ... S 8im), where (p) is the p-th new index after sorting. Then, the "Sparsity" of g is defined by
Sparsity(3) =
Σ || p=1
2
P
(3mml
P -- (P
(133)
P
This measure is based on weighted sum of all coefficients, evaluating how important particular coefficient is to an overall sparsity. We have Sparsity(3) = 0 for the least sparse case where all coefficients have the same value, and Sparsity(S) = 1 for the most sparse case where there exist few non-zero coefficients under the situation p co. Therefore, the measure can evaluate sparsity of g.
1.2 Comparison with Other Methods
Adding to Section I. we also compared the proposed model with a linear regression model without regular- ization, assuming the Azzalini's skew-normal distribution, skew-f distribution, and skew-Cauchy distribution for noise (Arellano-Vallee Valle & Azzalim, 2013; Azzalimi & Arellano-Valle 2013 Azzalini & Salehi, 2020). These methods are denoted by "skew-N," "skew-t," and "skew-C," " respectively. The experimental setting and data are the same as in Section [T.
Four evaluation measures, i.e., MSE()), MSE(3), Sparsity($), and Model Size($), are shown in Figure F. The definitions of these measures are the same as in Section A. As result of Figure B. the proposed method outperformed the comparative methods in all cases, indicating better prediction with a smaller number of features. Note that the model size for "skew-N," "skew-t," " and "skew-C" always 100 due to no regularization term, and the results of "Proposed," "Maeedd"""""C "Chen+," and "Lasso" are the same as in Figure A
1.3 Real-World Financial Data
We applied the proposed method the Engineering Graduate Salary (EGS) prediction data (Aggarwa atl MULE), which provides engineering graduates' employment outcomes with standardized assessment scores. EGS has fewer features than the two medical datasets in Section вak, and the meanings of all features are specifically given. These properties allow to consider the validity of the estimated active features. For simplicity, we used only numerical type features (sce Table E for details). The data used here consisted of 3998 and P 25. Other experimental conditions were the same as in Section B
Table shows the result EGS with the means and standard deviations (in parentheses) of the evaluation measures. The proposod method simultancously achieved the smallest. prediction error and the smallest model size. Regarding the residuals of the test data after training with Lasso, the mean of the normalized skewness of the residuals was 2.24.
We also obtained interesting results. Figure B shows box plots of the estimated coefficients for cach method. The top 10 coefficients are picked out based on the averaged absolute values over 30 trials. For example, a graduate's quantitative ability ("Quant") ranked first. for the proposed method, while a graduation year ("12graduation") ranked first for Lasso and Lasso+YJJ fourth with about half an effect. graduation year. As the box plots indicate, the coefficient. of this feature sometimes zero. Moreover, an English score ("English") and college tier ("College ranked second and third in the proposed method. From these points, the proposed method could select more reasonable features.
the proposed method, "12graduation" ranked may sound strange that the most relevant factor for salary is the
32

---

Published in Transactions on Machine Learning Research (09/2024)
MSE
MSE (8)
Sparsity
Model Size (8)
2.2
2.0
2.0
6.7
1.6
1.0
0.0
1.5
0.5
0.0
0.4
2.5
-100
1.0
1.00
2.0
0.5
3.4
0.8
1.5
0.7
3.2
1.0
0.6
3.0
X16
MSE(55)
1.0
Sparsity
1.00
(8)
0.9
0.8
6.0
5.5
5.0
MSE ()]
MSE (5)
Sparsity (8)
Model Size (5)
1.00
2.3
1.5
0.5
0.4
Propesed
Chenk
Lacso
Skrae
Scann
Proporeed
Coentt
Lucsa
Stms-N
Spemd
StamiC
Proposed
Coant
Loc3o
Skow!
Stowi
Proposen
Chane
Lacso
Snowe
Figure 5: Comparison of six methods with four different types of simulation noises in terms of four evaluation measures. From top to bottom, each results from Data 1, Data 2, Data 3, and Data 4, respectively. All experiments were conducted for 50 runs with different random seeds.
Table 3: Results of EGS (N = 3998, P = 25). MSE(ÿ)
Sparsity($))
Size (9)
Proposed
2.11 x 10a (3.14 x 10¹)
4.62 x 10 (2.84 10-
2.31 x 10a
(1.05 10")
Chen+
4.97 x 10a (7.59 x 109)
1.16 x 10
(16.8 x 10-2)
8.33 x 10a (12.0 x 10")
Lasso
4.89 X 10 (7.57 x 10')
3.78 10-1
(3.72 x 10
2.44
10a
(0.68 x 10%)
Lasso + YJ
2.64 x 10a (3.86 x 101)
3.78 10-1 (3.72 x
2.44 x 10a
(0.68 x 10")
33

---

Published in Transactions on Machine Learning Research (09/2024)
Proposed
Chen+
Quant: #1
GraduationYear: #1
English: #2
CollegeTier: #2
CollegeTier: #3
HHH+
English: #3
12graduation: #4
Quant: #4
10percentage: #5
+8+
12percentage:
#5
12percentage: #6
agreeableness:
#6
ComputerScience:
+IH
conscientiousness: #7
Domain:
+
computerrogramminn: #8
collegeGPA:
+[+
collegeGPA:
GraduationYear:
Lasso
Lasso+Y)
12graduation:
12graduation: #1
CollegeTier:
CollegeTier: #2
+I--0
12percentage:
1Zpercentage:
#3
HIH
English:
English: #4
Quant: #5
HI+0
Quant: #5
H1-
ComputerScience: #6
ComputerScience: #6
+E+
10percentage:
#7
10percentage:
#7
+0+
GraduationYear:
#8
GraduationYear:
#8
conscientiousnes::
#9
+0
conscientiousness:
#9
*|
agreeableness: #10
agreeableness:
#10
HBi
-0.2
-0.1
0.0
0.1
0.2
-0.2
-0.1
0.0
0.1
0.2
Figure 6: Box plots of the estimated coefficients for EGS. In the vertical axis, the top 10 coefficients are picked out based on the averaged absolute values over 30 runs with different random seeds. Details of each feature are described Table EL
Table 4: Description of features in EGS data (see Table in Aggarwal et al. (2010) for details). Feature
Description
10percentage
Overall marks obtained in grade 10 examinations Year of graduation (senior year high school)
2graduution
12percentage
Overall marks obtained in grade 12 examinations
CollegeID
Unique identifying the college which the candidate attended Tier of college (computed from the average AMCAT scores) Aggregate GPA at graduation A unique ID to identify the city in which the college is located The tier of the city in which the college is located Year of graduation (Bachelor's degree) Score in AMCAT's English section Score in AMCAT's Logical ability section
CollegeTiae Tier
CollegeGIA
College CityID
CollegeCity CityTicr Tier
Graduation Year
English
Logical
Quant
Score AMCAT's Quantitative ability section
Domain
Score in AMCAT's Domain module
ComputerrProraammig
Score in AMCAT's Computer programming section
Eleetronnes
Score in AMCAT's Electronics & Semiconduetor Engineering section Score in AMCAT's Computer Science section
ComputerScience
Mechanical Engg
Score in AMCAT's Mechanical Engineering section
ElectricalEnge
Score in AMCAT's Electrical Engineering section
TelecomEnge
Score AMCAT's Telecommmuication Engineering section
CivilEngg
Score in AMCAT's Civil Engineering section
conscientiousness
Score in one of the sections of AMCAT's personality test
agreeableness extraversion
Score in one of the sections of AMCAT's personality test
Score in one of the sections of AMCAT's personality test
mueroticism
Score in one of the sections of AMCAT's personality test
openess to experience
Score in one of the sections of AMCAT's personality test
34


Reviews:
[{"review_id": "FrrqFFTdv0", "reviewer": "Reviewer_j3Vk", "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\n- Novel and useful problem formulation. I believe the problem of regression with a skewed noise distribution is an important problem, and it is quite surprising that regression with a mode-invariant noise distribution has not been studied.\n- The authors’ proposed alternating minimization method is well-motivated.\n- Theoretical guarantees are interesting and non-trivial. In particular, I think Theorem 3 is a nice result.\n- Simulated results adequately demonstrate the improved performance of their method, particularly in the presence of model mis-specification.\n\nWeaknesses (see below for more comments):\n\n- Writing is somewhat unclear, e.g. definitions for variables are missing, language can be confusing.\n- Theory could be validated empirically.\n- Real data evaluation could be investigated more thoroughly.", "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": "In this manuscript the authors formulate a linear regression problem where the noise distribution is skewed and \"mode-invariant\", meaning that for any choice of the distribution parameters, the location parameter is equal to the mode. The authors present a novel alternating minimization method for solving this problem; one of the minimization problems is convex, while the other is non-convex but only involves two scalar variables. The authors prove theoretical guarantees on the excess risk, sparsity, and likelihood of their proposed solution, under certain assumptions. Finally, the authors demonstrate the improved performance of their method on simulated data -- including data simulated under model mis-specification -- and on two real datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Below I elaborate on the above weaknesses and what I would like to see.\n\n[Critical to securing recommendation.] Improving writing.\n\nThe writing of the manuscript can often be unclear, with definitions sometimes not given and non-standard terminology used. Below I list some specific examples that I would like to see updated.\n\n- Abstract: authors note that “regression model is always for a mode”. I have not heard this terminology before and it is imprecise. I suggest the authors change this to instead note that the noise distribution  always  has an interpretable mode parameter.\n    - e.g. similarly, in the introduction, the authors note that “Chen et al treated a regression model on…” which is also imprecise terminology.\n    - Sec 2.4: “that regresses the location parameter” also is imprecise.\n- Figure 1: label the blue/black lines in the caption.\n- Some variables are not defined and some equations are not clearly explained.\n    - $f_{\\psi}$ is not defined in (13)\n    - Sec 3.2: authors should explicitly state what $\\psi^0$ and $\\psi$ are. Is $\\psi^0$ the true model?\n    - Sec 3.3: “has the order slightly larger than…”. I am not sure what “has the order” means. I suspect the precise terminology you want is that the excess risk/sparsity of your method is (log N)^2 larger than that of ordinary Lasso?\n        - Same comment for Sec 3.4.\n    - Eq (58): $v_n$ is introduced without definition. I assume it is equal to $r_{\\alpha}(u_n)$?\n    - Moreover, eq (59) does not immediately follow from Prop 5. The authors should write out this argument more clearly.\n\n[Would strengthen manuscript]\n\n- Theory could be validated empirically.\n\nThe authors make several assumptions in their convergence proofs, e.g. the restricted eigenvalue condition, and also present theoretical results on the quadratic margin (Theorem 3) and convergence (Theorem 4). It would be useful to see these results validated empirically, e.g. how loose/tight are the bounds derived by the authors?\n\n- Real data evaluation is weak.\n\nWhile the real data evaluation is interesting, the authors only perform a shallow comparison of their method versus others by looking at summary statistics. I suggest the authors compare the results of different models in more detail. For example, the authors should compare which coefficients are found by the authors approach that are not found by others, and vice-versa. Are the coefficients uniquely found by the authors’ approach more medically plausible? Do Lasso/other approaches find many false positives?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "OVCK7AdG3x", "reviewer": "Reviewer_7Eky", "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": "S: the problem where the noise in the data is skewed is both interesting and practically relevant. The method allows finding the level of skewedness in the data, which is a plus\n\nW: I think the main weakness is the optimization proposal lacks motivation and needs additional comparison.\nW1. Motivation: the main motivation is stated as \"the solution could be unstable in optimization due to the non-convexity\n of the Lasso-type problem.\" This is not true. See https://arxiv.org/abs/2210.01212, where the authors proved that L1 regularization can be solved efficiently for generic nonlinear objective functions with simple gradient descent\n\nW2. This is related to W1. I believe that the authors should avoid reinventing the wheel. As generic nonlinear solvers of L1 exists, the authors should either (1) rely on existing optimization method or (2) propose a new method and show that it has advantage over existing methods. \n\nSo far, the manuscript only proposes an optimization method without any motivation for the need of it, and I think this point needs to be improved", "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 to minimize the log-probability of a skewed distribution in order to fit a sparse linear regression model with skewed noise", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See the two weakness. I believe the authors need to benchmark the proposed optimization method against known ones or remove any claim of contribution based on it", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "XqJ3N0nm8Q", "reviewer": "Reviewer_rUk6", "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.**\n1. The result is novel and complete. The noise model considered is different from that in the previous work by Chen et al. (2014). Both the computational and statistical aspects are discussed. \n\n**Weaknesses.**\n1. The paper argues that the considered noise model is more interpretable. However, interpretability is a vague notion and lacks an objective definition. It would be beneficial if the authors could provide more objective, convincing motivation for their work.\n2. The optimization method for computing the estimate follows a majorization-minimization argument. However, the method appears to be heuristic and lacks a convergence guarantee. \n3. Theorem 4, which ensures that the guarantees of the estimation error and variable screening ability hold with high probability, is obtained by applying the results of Bühlmann & Van De Geer (2011), which were derived for the FMR model. The FMR model is a mixture-of-Gaussian linear regression model. It is unclear to me why the linear regression model with skewed noise, which this paper studies, is considered a special case of the FMR model.\n4. Below are some minor comments. \n    - The function $H_\\alpha$ looks confusing in (3), as its definition appears far from (3) in the next paragraph. \n    - Last line on p. 4: The word \"bound\" should be replaced with \"bounded.\"\n    - The purpose of Section 3.1 is unclear from the main text.\n    - Second paragraph of Section 3.1: The meaning of the notation $\\mathcal{X} = \\{ X \\}$ is unclear to me. \n    - Assumption 2: Define $\\kappa = 1 / \\tilde{\\kappa}$ can be confusing.", "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 paper studies the $\\ell_1$-penalized maximum-likelihood estimator for sparse linear regression with skewed noise. Whereas the same topic has been studied by Chen et al. (2014), this paper adopts a different noise model called the mode-invariant skew-normal distribution, proposed by Fujisawa and Abe (2015). A heuristic optimization method to compute the $\\ell_1$-penalized maximum-likelihood (ML) estimate is provided. Both the estimation error and the variable screening ability of the estimator are analyzed.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please address the weaknesses above.", "broader_impact_concerns": "N/A.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper worked on a difficult regression problem with a linear regression model but the noise can be highly skewed (showed in Eq. (5)). It was well motivated to adopt a certain mode-invariant distribution as the model of the noise (showed in Fig. 1). There were also some theoretical guarantees for the proposed method. After the rebuttal, all the three reviewers were positive, and thus I think we can accept the paper for publication. Please incorporate your rebuttal into the camera-ready version.
Decision: Accept (poster)