11 #include <shogun/lib/config.h> 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 >
get_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);
316 map_beta = (1-s)*map_beta_prev + s*map_beta;
325 deactivate_variable(i_kick);
328 int32_t numRows = map_Xa.rows();
329 int32_t numCols = map_Xa.cols()-1;
330 if( i_kick < numCols )
331 map_Xa.block(0, i_kick, numRows, numCols-i_kick) =
332 map_Xa.block(0, i_kick+1, numRows, numCols-i_kick).eval();
336 m_beta_path_t.push_back(beta);
338 m_beta_idx.push_back(nloop);
340 m_beta_idx[m_num_active] = nloop;
345 SG_DEBUG(
"Added : %d , Dropped %d, Active set size %d max_corr %.17f \n", i_max_corr, i_kick, m_num_active, max_corr);
349 for(
size_t i = 0; i < m_beta_path_t.size(); ++i)
351 std::vector<float64_t> va;
352 for(
size_t p = 0; p < m_beta_path_t[i].size(); ++p){
353 va.push_back((
float64_t) m_beta_path_t[i][p]);
355 m_beta_path.push_back(va);
365 template <
typename ST>
371 ST diag_k = map_X.col(i_max_corr).dot(map_X.col(i_max_corr));
379 map_R.transpose().template triangularView<Lower>().
template solveInPlace<OnTheLeft>(R_k);
385 map_R_new.block(0, 0, num_active, num_active) = map_R;
386 sg_memcpy(R_new.matrix+num_active*(num_active+1), R_k.
data(),
sizeof(ST)*(num_active));
387 map_R_new.row(num_active).setZero();
388 map_R_new(num_active, num_active) = R_kk;
392 template <
typename ST>
395 if (i_kick != m_num_active-1)
398 for (
index_t j=i_kick; j < m_num_active-1; ++j)
399 for (
index_t i=0; i < m_num_active; ++i)
403 for (
index_t i=i_kick; i < m_num_active-1; ++i)
405 plane_rot(R(i,i),R(i+1,i), R(i,i), R(i+1,i), G);
406 if (i < m_num_active-2)
408 for (
index_t k=i+1; k < m_num_active-1; ++k)
411 ST Rik = R(i,k), Ri1k = R(i+1,k);
412 R(i,k) = G(0,0)*Rik + G(0,1)*Ri1k;
413 R(i+1,k) = G(1,0)*Rik+G(1,1)*Ri1k;
420 for (
index_t i=0; i < m_num_active-1; ++i)
421 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)
SGMatrix< ST > get_feature_matrix()
virtual ELabelType get_label_type() const =0
Real Labels are real-valued labels.
virtual void set_w(const SGVector< float64_t > src_w)
virtual int32_t get_num_labels() const =0
real valued labels (e.g. for regression, classifier outputs)
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)
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 float64_t onenorm(T *x, int32_t len)
|| x ||_1
static void plane_rot(ST x0, ST x1, ST &y0, ST &y1, SGMatrix< ST > &G)
void switch_w(int32_t num_variable)
T * get_column_vector(index_t col) const
int32_t get_max_non_zero() const
virtual EFeatureClass get_feature_class() const =0
Class LinearMachine is a generic interface for all kinds of linear machines like classifiers.
bool train_machine(CFeatures *data)
void set_epsilon(float64_t epsilon)
int32_t get_num_features() const
virtual int32_t get_num_vectors() const
all of classes and functions are contained in the shogun namespace
The class Features is the base class of all feature objects.
int32_t get_path_size() const
static float32_t sqrt(float32_t x)
T max(const Container< T > &a)
float64_t get_max_l1_norm() const
virtual EFeatureType get_feature_type() const =0