23 using namespace Eigen;
27 CLeastAngleRegression::CLeastAngleRegression(
bool lasso) :
29 m_max_nonz(0), m_max_l1_norm(0)
42 template <
typename ST>
44 int32_t &imax, ST& vmax)
48 for (
size_t i=0; i < vec.size(); ++i)
61 template <
typename ST>
69 G(0, 0) = G(1, 1) = 1;
117 SG_ERROR(
"Feature-type (%d) must be of type F_SHORTREAL (%d), F_DREAL (%d) or F_LONGREAL (%d).\n",
123 template <
typename ST>
126 std::vector<std::vector<ST>> m_beta_path_t;
131 bool lasso_cond =
false;
132 bool stop_cond =
false;
136 m_beta_path_t.clear();
138 m_active_set.clear();
139 m_is_active.resize(n_fea);
140 fill(m_is_active.begin(), m_is_active.end(),
false);
151 map_X = map_Xr.transpose();
156 vector<ST> beta(n_fea);
158 vector<ST> Xy(n_fea);
164 vector<ST> mu(n_vec);
167 vector<ST> corr(n_fea);
169 vector<ST> corr_sign(n_fea);
175 int32_t i_max_corr = 1;
178 m_beta_path_t.push_back(beta);
179 m_beta_idx.push_back(0);
182 int32_t max_active_allowed =
CMath::min(n_vec-1, n_fea);
188 while (m_num_active < max_active_allowed && max_corr/n_vec > m_epsilon && !stop_cond)
194 map_corr = map_Xy - (map_Xr*map_mu);
197 for (
size_t i=0; i < corr.size(); ++i)
206 if (m_num_active == 0)
210 ST diag_k = map_X.col(i_max_corr).dot(map_X.col(i_max_corr));
215 activate_variable(i_max_corr);
221 map_Xa.col(m_num_active-1)=map_X.col(i_max_corr);
224 for (
index_t i=0; i < m_num_active; ++i)
225 corr_sign_a[i] = corr_sign[m_active_set[i]];
229 typename SGVector<ST>::EigenVectorXt solve = map_R.transpose().template triangularView<Lower>().
template solve<OnTheLeft>(map_corr_sign_a);
234 ST AA = GA1.dot(map_corr_sign_a);
245 ST gamma = max_corr / AA;
246 if (m_num_active < n_fea)
248 #pragma omp parallel for shared(gamma)
249 for (
index_t i=0; i < n_fea; ++i)
255 ST dir_corr = map_u.dot(map_X.col(i));
257 ST tmp1 = (max_corr-corr[i])/(AA-dir_corr);
258 ST tmp2 = (max_corr+corr[i])/(AA+dir_corr);
270 int32_t i_change=i_max_corr;
277 for (
index_t i=0; i < m_num_active; ++i)
279 ST tmp = -beta[m_active_set[i]] / wA(i);
287 if (lasso_bound < gamma)
291 i_change = m_active_set[i_kick];
296 map_mu += gamma*map_u;
299 for (
index_t i=0; i < m_num_active; ++i)
300 beta[m_active_set[i]] += gamma * wA(i);
302 if (m_max_l1_norm > 0)
305 if (l1 > m_max_l1_norm)
311 ST s = (m_max_l1_norm-l1_prev)/(l1-l1_prev);
315 map_beta = (1-s)*map_beta_prev + s*map_beta;
324 deactivate_variable(i_kick);
327 int32_t numRows = map_Xa.rows();
328 int32_t numCols = map_Xa.cols()-1;
329 if( i_kick < numCols )
330 map_Xa.block(0, i_kick, numRows, numCols-i_kick) =
331 map_Xa.block(0, i_kick+1, numRows, numCols-i_kick).eval();
335 m_beta_path_t.push_back(beta);
336 if (
size_t(m_num_active) >= m_beta_idx.size())
337 m_beta_idx.push_back(nloop);
339 m_beta_idx[m_num_active] = nloop;
342 if (m_max_nonz > 0 && m_num_active >= m_max_nonz)
344 SG_DEBUG(
"Added : %d , Dropped %d, Active set size %d max_corr %.17f \n", i_max_corr, i_kick, m_num_active, max_corr);
348 for(
size_t i = 0; i < m_beta_path_t.size(); ++i)
350 std::vector<float64_t> va;
351 for(
size_t p = 0; p < m_beta_path_t[i].size(); ++p){
352 va.push_back((
float64_t) m_beta_path_t[i][p]);
354 m_beta_path.push_back(va);
364 template <
typename ST>
370 ST diag_k = map_X.col(i_max_corr).dot(map_X.col(i_max_corr));
378 map_R.transpose().template triangularView<Lower>().
template solveInPlace<OnTheLeft>(R_k);
384 map_R_new.block(0, 0, num_active, num_active) = map_R;
385 memcpy(R_new.matrix+num_active*(num_active+1), R_k.
data(),
sizeof(ST)*(num_active));
386 map_R_new.row(num_active).setZero();
387 map_R_new(num_active, num_active) = R_kk;
391 template <
typename ST>
394 if (i_kick != m_num_active-1)
397 for (
index_t j=i_kick; j < m_num_active-1; ++j)
398 for (
index_t i=0; i < m_num_active; ++i)
402 for (
index_t i=i_kick; i < m_num_active-1; ++i)
404 plane_rot(R(i,i),R(i+1,i), R(i,i), R(i+1,i), G);
405 if (i < m_num_active-2)
407 for (
index_t k=i+1; k < m_num_active-1; ++k)
410 ST Rik = R(i,k), Ri1k = R(i+1,k);
411 R(i,k) = G(0,0)*Rik + G(0,1)*Ri1k;
412 R(i+1,k) = G(1,0)*Rik+G(1,1)*Ri1k;
419 for (
index_t i=0; i < m_num_active-1; ++i)
420 for (
index_t j=0; j < m_num_active-1; ++j)
virtual const char * get_name() const =0
static const float64_t MACHINE_EPSILON
SGMatrix< ST > cholesky_delete(SGMatrix< ST > &R, int32_t i_kick)
virtual ELabelType get_label_type() const =0
Real Labels are real-valued labels.
int32_t get_num_features() const
virtual int32_t get_num_labels() const =0
real valued labels (e.g. for regression, classifier outputs)
SGMatrix< ST > get_feature_matrix()
virtual int32_t get_num_vectors() const =0
virtual ~CLeastAngleRegression()
static void find_max_abs(const std::vector< ST > &vec, const std::vector< bool > &ignore_mask, int32_t &imax, ST &vmax)
static void plane_rot(ST x0, ST x1, ST &y0, ST &y1, SGMatrix< ST > &G)
void switch_w(int32_t num_variable)
virtual int32_t get_num_vectors() const
virtual EFeatureClass get_feature_class() const =0
T * get_column_vector(index_t col) const
Class LinearMachine is a generic interface for all kinds of linear machines like classifiers.
bool train_machine(CFeatures *data)
all of classes and functions are contained in the shogun namespace
SGMatrix< ST > cholesky_insert(const SGMatrix< ST > &X, const SGMatrix< ST > &X_active, SGMatrix< ST > &R, int32_t i_max_corr, int32_t num_active)
The class Features is the base class of all feature objects.
static float64_t onenorm(T *x, int32_t len)
|| x ||_1
Matrix::Scalar max(Matrix m)
static float32_t sqrt(float32_t x)
virtual EFeatureType get_feature_type() const =0