25 using namespace shogun;
30 vector<float64_t> result(size);
31 fill(result.begin(), result.end(), val);
41 G(0, 0) = G(1, 1) = 1;
61 static void find_max_abs(
const vector<float64_t> &vec,
const vector<bool> &ignore,
66 for (uint32_t i=0; i < vec.size(); ++i)
80 m_max_nonz(0), m_max_l1_norm(0)
96 SG_ERROR(
"Expected RegressionLabels\n");
105 SG_ERROR(
"Number of training vectors does not match number of labels\n");
108 SG_ERROR(
"Expected Simple Features\n");
111 SG_ERROR(
"Expected Real Features\n");
117 bool lasso_cond =
false;
118 bool stop_cond =
false;
124 m_active_set.clear();
125 m_is_active.resize(n_fea);
126 fill(m_is_active.begin(), m_is_active.end(),
false);
135 for (int32_t i=0; i < n_vec; ++i)
137 for (int32_t j=0; j < n_fea; ++j)
146 cblas_dgemv(CblasColMajor, CblasTrans, n_vec, n_fea, 1, X.
matrix, n_vec,
147 y.vector, 1, 0, &Xy[0], 1);
155 vector<float64_t> corr_sign(n_fea);
161 int32_t i_max_corr = 1;
164 m_beta_path.push_back(beta);
165 m_beta_idx.push_back(0);
174 copy(Xy.begin(), Xy.end(), corr.begin());
175 cblas_dgemv(CblasColMajor, CblasTrans, n_vec, n_fea, -1,
176 X.
matrix, n_vec, &mu[0], 1, 1, &corr[0], 1);
179 for (uint32_t i=0; i < corr.size(); ++i)
188 R=cholesky_insert(X, R, i_max_corr);
189 activate_variable(i_max_corr);
193 vector<float64_t> corr_sign_a(m_num_active);
194 for (int32_t i=0; i < m_num_active; ++i)
195 corr_sign_a[i] = corr_sign[m_active_set[i]];
198 vector<float64_t> GA1(corr_sign_a);
199 cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit,
200 m_num_active, 1, 1, R.
matrix, m_num_active, &GA1[0], m_num_active);
201 cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit,
202 m_num_active, 1, 1, R.
matrix, m_num_active, &GA1[0], m_num_active);
205 float64_t AA = cblas_ddot(m_num_active, &GA1[0], 1, &corr_sign_a[0], 1);
209 vector<float64_t> wA(GA1);
210 for (int32_t i=0; i < m_num_active; ++i)
216 for (int32_t i=0; i < m_num_active; ++i)
219 cblas_daxpy(n_vec, wA[i],
225 if (m_num_active < n_fea)
227 for (int32_t i=0; i < n_fea; ++i)
236 float64_t tmp1 = (max_corr-corr[i])/(AA-dir_corr);
237 float64_t tmp2 = (max_corr+corr[i])/(AA+dir_corr);
246 int32_t i_change=i_max_corr;
251 float64_t lasso_bound = numeric_limits<float64_t>::max();
253 for (int32_t i=0; i < m_num_active; ++i)
255 float64_t tmp = -beta[m_active_set[i]] / wA[i];
263 if (lasso_bound < gamma)
267 i_change = m_active_set[i_kick];
272 cblas_daxpy(n_vec, gamma, &u[0], 1, &mu[0], 1);
275 for (int32_t i=0; i < m_num_active; ++i)
276 beta[m_active_set[i]] += gamma * wA[i];
279 if (m_max_l1_norm > 0)
282 if (l1 > m_max_l1_norm)
288 float64_t s = (m_max_l1_norm-l1_prev)/(l1-l1_prev);
292 cblas_dscal(n_fea, s, &beta[0], 1);
293 cblas_daxpy(n_fea, 1-s, &m_beta_path[nloop][0], 1, &beta[0], 1);
303 R=cholesky_delete(R, i_kick);
304 deactivate_variable(i_kick);
308 m_beta_path.push_back(beta);
309 if (
size_t(m_num_active) >= m_beta_idx.size())
310 m_beta_idx.push_back(nloop);
312 m_beta_idx[m_num_active] = nloop;
315 if (m_max_nonz > 0 && m_num_active >= m_max_nonz)
334 if (m_num_active == 0)
344 vector<float64_t> col_k(m_num_active);
345 for (int32_t i=0; i < m_num_active; ++i)
353 vector<float64_t> R_k(col_k);
354 cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, m_num_active, 1,
355 1, R.
matrix, m_num_active, &R_k[0], m_num_active);
358 cblas_ddot(m_num_active, &R_k[0], 1, &R_k[0], 1));
362 for (int32_t i=0; i < m_num_active; ++i)
363 for (int32_t j=0; j < m_num_active; ++j)
365 for (int32_t i=0; i < m_num_active; ++i)
366 nR(i, m_num_active) = R_k[i];
367 for (int32_t i=0; i < m_num_active; ++i)
368 nR(m_num_active, i) = 0;
369 nR(m_num_active, m_num_active) = R_kk;
378 if (i_kick != m_num_active-1)
381 for (int32_t j=i_kick; j < m_num_active-1; ++j)
382 for (int32_t i=0; i < m_num_active; ++i)
386 for (int32_t i=i_kick; i < m_num_active-1; ++i)
388 plane_rot(R(i,i),R(i+1,i), R(i,i), R(i+1,i), G);
389 if (i < m_num_active-2)
391 for (int32_t k=i+1; k < m_num_active-1; ++k)
395 R(i,k) = G(0,0)*Rik + G(0,1)*Ri1k;
396 R(i+1,k) = G(1,0)*Rik+G(1,1)*Ri1k;
403 for (int32_t i=0; i < m_num_active-1; ++i)
404 for (int32_t j=0; j < m_num_active-1; ++j)
410 #endif // HAVE_LAPACK