Skip to content

Commit 377a8f9

Browse files
statmlbenstatmlben
authored
Merge pull request #47 from DataboyUsen/main
Adaptive lasso (#10)
2 parents 7545e3c + 065450d commit 377a8f9

8 files changed

Lines changed: 402 additions & 52 deletions

File tree

‎doc/source/example.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,3 +32,4 @@ Example Gallery
3232
examples/GridSearchCV_reg_losses.ipynb
3333

3434
examples/ElasticNet.ipynb
35+
examples/Adaptive_ElasticNet.ipynb

‎doc/source/examples/Adaptive_ElasticNet.ipynb‎

Lines changed: 201 additions & 0 deletions
Large diffs are not rendered by default.

‎rehline/_base.py‎

Lines changed: 17 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -195,7 +195,7 @@ def call_ReLHLoss(self, score):
195195
if self.H > 0:
196196
rehu_input = (self._S.T * score[:, np.newaxis]).T + self._T
197197
return np.sum(_relu(relu_input), 0) + np.sum(_rehu(rehu_input, self._Tau), 0)
198-
198+
199199
@abstractmethod
200200
def fit(self, X, y, sample_weight):
201201
"""Fit model."""
@@ -286,7 +286,7 @@ def ReHLine_solver(
286286
T=None,
287287
A=None,
288288
b=None,
289-
rho=0.0,
289+
rho=None,
290290
Lambda=None,
291291
Gamma=None,
292292
xi=None,
@@ -307,6 +307,8 @@ def ReHLine_solver(
307307
A = np.empty(shape=(0, 0))
308308
if b is None:
309309
b = np.empty(shape=(0))
310+
if rho is None:
311+
rho = np.empty(shape=(0))
310312
if Lambda is None:
311313
Lambda = np.empty(shape=(0, 0))
312314
if Gamma is None:
@@ -330,14 +332,14 @@ def ReHLine_solver(
330332
X,
331333
A,
332334
b,
335+
rho,
333336
U,
334337
V,
335338
S,
336339
T,
337340
Tau,
338341
max_iter,
339342
tol,
340-
rho,
341343
shrink,
342344
verbose,
343345
trace_freq,
@@ -773,32 +775,32 @@ def _cast_sample_weight(U, V, Tau, S, T, C=1.0, sample_weight=None):
773775
# U_new = np.zeros((self.L+2, n+d))
774776
# V_new = np.zeros((self.L+2, n+d))
775777
# ## Block 1
776-
# if len(self.U):
777-
# U_new[:self.L, :n] = self.U
778-
# V_new[:self.L, :n] = self.V
778+
# if len(self._U):
779+
# U_new[:self.L, :n] = self._U
780+
# V_new[:self.L, :n] = self._V
779781
# ## Block 2
780782
# U_new[-2,n:] = l1_pen
781783
# U_new[-1,n:] = -l1_pen
782784

783-
# if len(self.S):
785+
# if len(self._S):
784786
# S_new = np.zeros((self.H, n+d))
785787
# T_new = np.zeros((self.H, n+d))
786788
# Tau_new = np.zeros((self.H, n+d))
787789

788-
# S_new[:,:n] = self.S
789-
# T_new[:,:n] = self.T
790-
# Tau_new[:,:n] = self.Tau
790+
# S_new[:,:n] = self._S
791+
# T_new[:,:n] = self._T
792+
# Tau_new[:,:n] = self._Tau
791793

792-
# self.S = S_new
793-
# self.T = T_new
794-
# self.Tau = Tau_new
794+
# self._S = S_new
795+
# self._T = T_new
796+
# self._Tau = Tau_new
795797

796798
# ## fake X
797799
# X_fake = np.zeros((n+d, d))
798800
# X_fake[:n,:] = X
799801
# X_fake[n:,:] = np.identity(d)
800802

801-
# self.U = U_new
802-
# self.V = V_new
803+
# self._U = U_new
804+
# self._V = V_new
803805
# self.auto_shape()
804806
# return X_fake

‎rehline/_class.py‎

Lines changed: 25 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -493,7 +493,7 @@ class plqERM_ElasticNet(_BaseReHLine, BaseEstimator):
493493
494494
.. math::
495495
496-
\min_{\mathbf{\beta} \in \mathbb{R}^d} C \sum_{i=1}^n \text{PLQ}(y_i, \mathbf{x}_i^T \mathbf{\beta}) + \text{l1_ratio} \| \mathbf{\beta} \|_1 + \frac{1}{2} (1 - \text{l1_ratio}) \| \mathbf{\beta} \|_2^2, \ \text{ s.t. } \
496+
\min_{\mathbf{\beta} \in \mathbb{R}^d} C \sum_{i=1}^n \text{PLQ}(y_i, \mathbf{x}_i^T \mathbf{\beta}) + \text{l1_ratio} \sum_{j=1}^d \omega_j | \beta_j | + \frac{1}{2} (1 - \text{l1_ratio}) \| \mathbf{\beta} \|_2^2, \ \text{ s.t. } \
497497
\mathbf{A} \mathbf{\beta} + \mathbf{b} \geq \mathbf{0},
498498
499499
The function supports various loss functions, including:
@@ -527,6 +527,9 @@ class plqERM_ElasticNet(_BaseReHLine, BaseEstimator):
527527
The ElasticNet mixing parameter, with 0 <= l1_ratio < 1. For l1_ratio = 0 the penalty
528528
is an L2 penalty. For 0 < l1_ratio < 1, the penalty is a combination of L1 and L2.
529529
530+
omega : array of shape (n_features, ), default=np.empty(shape=(0, 0))
531+
Weight coefficients for adaptive lasso.
532+
530533
verbose : int, default=0
531534
Enable verbose output. Note that this setting takes advantage of a
532535
per-process runtime setting in liblinear that, if enabled, may not work
@@ -584,6 +587,7 @@ def __init__(
584587
constraint=None,
585588
C=1.0,
586589
l1_ratio=0.5,
590+
omega=None,
587591
U=None,
588592
V=None,
589593
Tau=None,
@@ -602,8 +606,8 @@ def __init__(
602606
self.constraint = constraint if constraint is not None else []
603607
self.C = C
604608
self.l1_ratio = l1_ratio
605-
self.rho = l1_ratio / (1 - l1_ratio)
606609
self.C_eff = C / (1 - l1_ratio)
610+
self.omega = omega if omega is not None else np.empty(shape=(0, 0))
607611
self._U = U if U is not None else np.empty(shape=(0, 0))
608612
self._V = V if V is not None else np.empty(shape=(0, 0))
609613
self._S = S if S is not None else np.empty(shape=(0, 0))
@@ -678,6 +682,25 @@ def fit(self, X, y, sample_weight=None):
678682
self._xi = np.empty(shape=(0, 0))
679683
self._mu = np.empty(shape=(0, 0))
680684

685+
if self.l1_ratio == 0:
686+
self.rho = None
687+
if self.omega.size > 0:
688+
warnings.warn(
689+
f"Omega will be ignored since l1_ratio=0.",
690+
UserWarning,
691+
stacklevel=2
692+
)
693+
else:
694+
if self.omega.size not in (0, d):
695+
raise ValueError(
696+
f"Omega length {self.omega.size} must be 0 or {d} (n_features)"
697+
)
698+
if not np.all(self.omega > 0):
699+
raise ValueError(
700+
"All elements in omega must be strictly positive."
701+
)
702+
self.rho = np.full(d, self.l1_ratio / (1 - self.l1_ratio)) * (self.omega if self.omega.size == d else 1.0)
703+
681704
result = ReHLine_solver(
682705
X=X,
683706
U=U_weight,

‎rehline/_internal.pyi‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,14 +19,14 @@ def rehline_internal(
1919
X: npt.NDArray[np.float64],
2020
A: npt.NDArray[np.float64],
2121
b: npt.NDArray[np.float64],
22+
rho: npt.NDArray[np.float64],
2223
U: npt.NDArray[np.float64],
2324
V: npt.NDArray[np.float64],
2425
S: npt.NDArray[np.float64],
2526
T: npt.NDArray[np.float64],
2627
Tau: npt.NDArray[np.float64],
2728
max_iter: int,
2829
tol: float,
29-
rho: float = ...,
3030
shrink: int = ...,
3131
verbose: int = ...,
3232
trace_freq: int = ...,

‎src/rehline.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -19,15 +19,15 @@ using ReHLineResult = rehline::ReHLineResult<Matrix>;
1919

2020
void rehline_internal(
2121
ReHLineResult& result,
22-
const MapMat& X, const MapMat& A, const MapVec& b,
22+
const MapMat& X, const MapMat& A, const MapVec& b, const MapVec& rho,
2323
const MapMat& U, const MapMat& V,
2424
const MapMat& S, const MapMat& T, const MapMat& Tau,
25-
int max_iter, double tol, double rho = 0.0, int shrink = 1,
25+
int max_iter, double tol, int shrink = 1,
2626
int verbose = 0, int trace_freq = 100
2727
)
2828
{
29-
rehline::rehline_solver(result, X, A, b, U, V, S, T, Tau,
30-
max_iter, tol, rho, shrink, verbose, trace_freq);
29+
rehline::rehline_solver(result, X, A, b, rho, U, V, S, T, Tau,
30+
max_iter, tol, shrink, verbose, trace_freq);
3131
}
3232

3333
PYBIND11_MODULE(_internal, m) {

‎src/rehline.h‎

Lines changed: 34 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,7 @@ void reset_fv_set(std::vector<std::pair<Index, Index>>& fvset, std::size_t n, st
8080
// * S, T, Tau: [H x n]
8181
// * A : [K x d]
8282
// * b : [K]
83+
// * rho : [d]
8384
// - Pre-computed
8485
// * r: [n]
8586
// * p: [K]
@@ -139,6 +140,7 @@ class ReHLineSolver
139140
const Index m_L;
140141
const Index m_H;
141142
const Index m_K;
143+
const Index m_W;
142144

143145
// Input matrices and vectors
144146
RMatrix m_X;
@@ -149,7 +151,7 @@ class ReHLineSolver
149151
ConstRefMat m_Tau;
150152
RMatrix m_A;
151153
ConstRefVec m_b;
152-
Scalar m_rho; // l1_ratio / (1 - l1_ratio)
154+
ConstRefVec m_rho;
153155

154156
// Pre-computed
155157
Vector m_gk_denom; // ||a[k]||^2
@@ -174,7 +176,7 @@ class ReHLineSolver
174176
// =================== Initialization functions =================== //
175177

176178
// Compute the primal variable beta from dual variables
177-
// beta = A'xi - U3 * vec(Lambda) - S3 * vec(Gamma)
179+
// beta = A'xi - U3 * vec(Lambda) - S3 * vec(Gamma) + 2 * Mu - rho
178180
// A can be empty, one of U and V may be empty
179181
inline void set_primal()
180182
{
@@ -193,9 +195,9 @@ class ReHLineSolver
193195
if (m_H > 0)
194196
LHterm.noalias() += m_S.cwiseProduct(m_Gamma).colwise().sum().transpose();
195197
m_beta.noalias() -= m_X.transpose() * LHterm;
196-
197-
if (m_rho > 0)
198-
m_beta.noalias() += Scalar(2.0) * m_mu - m_rho * Vector::Ones(m_d);
198+
// L1 related
199+
if (m_W > 0)
200+
m_beta.noalias() += Scalar(2.0) * m_mu - m_rho;
199201
}
200202

201203
// =================== Evaluating objection function =================== //
@@ -224,7 +226,7 @@ class ReHLineSolver
224226
// Quadratic term
225227
result += Scalar(0.5) * m_beta.squaredNorm();
226228
// L1 penalty term
227-
result += m_rho * m_beta.template lpNorm<1>();
229+
result += m_beta.cwiseAbs().cwiseProduct(m_rho).sum();
228230
return result;
229231
}
230232

@@ -251,8 +253,8 @@ class ReHLineSolver
251253
}
252254
// 2 * Mu - rho, [d x 1]
253255
Vector MuR = Vector::Zero(m_d);
254-
if (m_rho > 0)
255-
MuR = Scalar(2.0) * m_mu - m_rho * Vector::Ones(m_d);
256+
if (m_W > 0)
257+
MuR = Scalar(2.0) * m_mu - m_rho;
256258
// Compute dual objective function value
257259
Scalar obj = Scalar(0);
258260
// If K = 0, all terms that depend on A, xi, or b will be zero
@@ -261,30 +263,31 @@ class ReHLineSolver
261263
// 0.5 * ||Atxi||^2 - Atxi' * U3L - Atxi' * S3G + Atxi' MuR + xi' * b
262264
const Scalar Atxi_U3L = (m_L > 0) ? (Atxi.dot(U3L)) : Scalar(0);
263265
const Scalar Atxi_S3G = (m_H > 0) ? (Atxi.dot(S3G)) : Scalar(0);
264-
const Scalar Atxi_MuR = (m_rho > 0) ? (Atxi.dot(MuR)) : Scalar(0);
266+
const Scalar Atxi_MuR = (m_W > 0) ? (Atxi.dot(MuR)) : Scalar(0);
265267
obj += Scalar(0.5) * Atxi.squaredNorm() - Atxi_U3L - Atxi_S3G + Atxi_MuR + m_xi.dot(m_b);
266268
}
267269
// If L = 0, all terms that depend on U, V, or Lambda will be zero
268270
if (m_L > 0)
269271
{
270272
// 0.5 * ||U3L||^2 + U3L' * S3G - U3L' * MuR - tr(Lambda * V')
271273
const Scalar U3L_S3G = (m_H > 0) ? (U3L.dot(S3G)) : Scalar(0);
272-
const Scalar U3L_MuR = (m_rho > 0) ? (U3L.dot(MuR)) : Scalar(0);
274+
const Scalar U3L_MuR = (m_W > 0) ? (U3L.dot(MuR)) : Scalar(0);
273275
obj += Scalar(0.5) * U3L.squaredNorm() + U3L_S3G - U3L_MuR -
274276
m_Lambda.cwiseProduct(m_V).sum();
275277
}
276278
// If H = 0, all terms that depend on S, T, or Gamma will be zero
277279
if (m_H > 0)
278280
{
279281
// 0.5 * ||S3G||^2 - S3G' * MuR + 0.5 * ||Gamma||^2 - tr(Gamma * T')
280-
const Scalar S3G_MuR = (m_rho > 0) ? (S3G.dot(MuR)) : Scalar(0);
282+
const Scalar S3G_MuR = (m_W > 0) ? (S3G.dot(MuR)) : Scalar(0);
281283
obj += Scalar(0.5) * S3G.squaredNorm() - S3G_MuR + Scalar(0.5) * m_Gamma.squaredNorm() -
282284
m_Gamma.cwiseProduct(m_T).sum();
283285
}
284-
// If rho = 0, all terms that depend on rho, or Mu will be zero
285-
if (m_rho > 0)
286-
obj += Scalar(2.0) * m_mu.squaredNorm() - Scalar(2.0) * m_rho * m_mu.sum() +
287-
Scalar(0.5) * m_d * m_rho * m_rho;
286+
// If W = 0, all terms that depend on rho or Mu will be zero
287+
if (m_W > 0)
288+
// 2.0 * ||Mu||^2 - 2.0 * rho' * Mu + 0.5 * ||rho||^2
289+
obj += Scalar(2.0) * m_mu.squaredNorm() - Scalar(2.0) * m_rho.dot(m_mu) +
290+
Scalar(0.5) * m_rho.squaredNorm();
288291
return obj;
289292
}
290293

@@ -368,7 +371,7 @@ class ReHLineSolver
368371
// Update mu and beta
369372
inline void update_mu_beta()
370373
{
371-
if (m_rho <= 0)
374+
if (m_W <= 0)
372375
return;
373376

374377
// Save original Mu
@@ -598,7 +601,7 @@ class ReHLineSolver
598601
// Overloaded version based on free variable set
599602
inline void update_mu_beta(std::vector<Index>& fv_set, Scalar& min_pg, Scalar& max_pg)
600603
{
601-
if (m_rho <= 0)
604+
if (m_W <= 0)
602605
return;
603606

604607
// Permutation
@@ -618,12 +621,13 @@ class ReHLineSolver
618621
for (auto j: fv_set)
619622
{
620623
const Scalar mu_j = m_mu[j];
624+
const Scalar rho_j = m_rho[j];
621625

622626
// Compute g_j
623627
const Scalar g_j = m_beta[j];
624628
// PG and shrink
625629
Scalar pg;
626-
const bool shrink = pg_mu(mu_j, g_j, m_rho, lb, ub, pg);
630+
const bool shrink = pg_mu(mu_j, g_j, rho_j, lb, ub, pg);
627631
if (shrink)
628632
continue;
629633

@@ -632,7 +636,7 @@ class ReHLineSolver
632636
min_pg = std::min(min_pg, pg);
633637
// Compute new mu_j
634638
const Scalar candid = mu_j - g_j * Scalar(0.5);
635-
const Scalar newmu = std::max(Scalar(0), std::min(m_rho, candid));
639+
const Scalar newmu = std::max(Scalar(0), std::min(rho_j, candid));
636640
// Update mu and beta
637641
m_mu[j] = newmu;
638642
m_beta[j] += Scalar(2.0) * (newmu - mu_j);
@@ -647,8 +651,9 @@ class ReHLineSolver
647651
ReHLineSolver(ConstRefMat X, ConstRefMat U, ConstRefMat V,
648652
ConstRefMat S, ConstRefMat T, ConstRefMat Tau,
649653
ConstRefMat A, ConstRefVec b,
650-
Scalar rho) :
651-
m_n(X.rows()), m_d(X.cols()), m_L(U.rows()), m_H(S.rows()), m_K(A.rows()),
654+
ConstRefVec rho) :
655+
m_n(X.rows()), m_d(X.cols()), m_L(U.rows()), m_H(S.rows()), m_K(A.rows()),
656+
m_W(rho.rows()), // check if l1 penalty is implemented
652657
m_X(X), m_U(U), m_V(V), m_S(S), m_T(T), m_Tau(Tau), m_A(A), m_b(b),
653658
m_rho(rho),
654659
m_gk_denom(m_K), m_gli_denom(m_L, m_n), m_ghi_denom(m_H, m_n),
@@ -691,10 +696,10 @@ class ReHLineSolver
691696
// Gamma.fill(std::min(1.0, 0.5 * tau));
692697
}
693698

694-
// Each element of Mu satisfies 0 <= mu_j <= rho,
699+
// Each element of Mu satisfies 0 <= mu_j <= rho_j,
695700
// and we use min(0.5 * rho, 1) to initialize (rho can be Inf)
696-
if (m_rho > 0)
697-
m_mu.setConstant( std::min(Scalar(0.5) * m_rho, Scalar(1)) );
701+
if (m_W > 0)
702+
m_mu.noalias() = (m_rho * Scalar(0.5)).cwiseMin(Scalar(1));
698703

699704
// Set primal variable based on duals
700705
set_primal();
@@ -746,15 +751,15 @@ class ReHLineSolver
746751
}
747752

748753

749-
if (m_rho > 0)
754+
if (m_W > 0)
750755
{
751756
// Check shape of warmstart parameters
752757
if (mu_ws.size() != m_d){
753758
throw std::invalid_argument("mu_ws must have size d");
754759
}
755760
// Check values of warmstart parameters
756-
if ((mu_ws.array() < Scalar(0)).any() || (mu_ws.array() > m_rho).any()){
757-
throw std::invalid_argument("mu_ws must be in [0, rho]");
761+
if ((mu_ws.array() < Scalar(0)).any() || (mu_ws.array() > m_rho.array()).any()){
762+
throw std::invalid_argument("mu_ws_j must be in [0, rho_j]");
758763
}
759764
m_mu = mu_ws;
760765
}
@@ -928,10 +933,10 @@ template <typename DerivedMat, typename DerivedVec, typename Index = int>
928933
void rehline_solver(
929934
ReHLineResult<typename DerivedMat::PlainObject, Index>& result,
930935
const Eigen::MatrixBase<DerivedMat>& X, const Eigen::MatrixBase<DerivedMat>& A,
931-
const Eigen::MatrixBase<DerivedVec>& b,
936+
const Eigen::MatrixBase<DerivedVec>& b, const Eigen::MatrixBase<DerivedVec>& rho,
932937
const Eigen::MatrixBase<DerivedMat>& U, const Eigen::MatrixBase<DerivedMat>& V,
933938
const Eigen::MatrixBase<DerivedMat>& S, const Eigen::MatrixBase<DerivedMat>& T, const Eigen::MatrixBase<DerivedMat>& Tau,
934-
Index max_iter, typename DerivedMat::Scalar tol, typename DerivedMat::Scalar rho = 0,
939+
Index max_iter, typename DerivedMat::Scalar tol,
935940
Index shrink = 1, Index verbose = 0, Index trace_freq = 100,
936941
std::ostream& cout = std::cout
937942
)

0 commit comments

Comments
 (0)