DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
decomp.h
Go to the documentation of this file.
1 
30 #ifndef DECOMP_H
31 #define DECOMP_H
32 
33 #include <utils.h>
34 #include <linalg.h>
35 
36 namespace spams
37 {
38 
39 /* **************************
40  * Greedy Forward Selection
41  * **************************/
42 
53 
54 template <typename T>
55 void omp(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
56  const int *L, const T* eps, const T* lambda, const bool vecL = false,
57  const bool vecEps = false, const bool Lambda=false, const int numThreads=-1,
58  Matrix<T>* path = NULL);
59 
60 template <typename T>
61 void omp_mask(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha, const Matrix<bool>& mask,
62  const int *L, const T* eps, const T* lambda, const bool vecL = false,
63  const bool vecEps = false, const bool Lambda=false, const int numThreads=-1,
64  Matrix<T>* path = NULL);
65 
67 template <typename T>
68 void coreORMP(Vector<T>& scores, Vector<T>& norm, Vector<T>& tmp,
69  Matrix<T>& Un, Matrix<T>& Undn, Matrix<T>& Unds, Matrix<T>& Gs,
70  Vector<T>& Rdn, const AbstractMatrix<T>& G, Vector<int>& ind,
71  Vector<T>& RUn, T& normX, const T* eps, const int* L, const T* lambda,
72  T* path = NULL);
73 
74 
76 template <typename T>
77 void coreORMPB(Vector<T>& RtD, const AbstractMatrix<T>& G, Vector<int>& ind,
78  Vector<T>& coeffs, T& normX, const int L, const T eps, const T lambda = 0);
79 
80 /* **************
81  * LARS - Lasso
82  * **************/
83 
88 enum constraint_type { L1COEFFS, L2ERROR, PENALTY, SPARSITY, L2ERROR2, PENALTY2};
89 
104 template <typename T>
105 void lasso(const Matrix<T>& X, const Matrix<T>& D,
106  SpMatrix<T>& spalpha,
107  int L, const T constraint, const T lambda2 = 0, constraint_type mode = PENALTY,
108  const bool pos = false, const bool ols = false, const int numThreads=-1,
109  Matrix<T>* path = NULL, const int length_path=-1);
110 
111 template <typename T>
112 void lasso(const Data<T>& X, const AbstractMatrix<T>& G, const AbstractMatrix<T>& DtX,
113  SpMatrix<T>& spalpha,
114  int L, const T constraint, constraint_type mode = PENALTY,
115  const bool pos = false, const bool ols = false, const int numThreads=-1,
116  Matrix<T>* path = NULL, const int length_path=-1);
117 
119 template <typename T>
120 void lasso2(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
121  int L, const T constraint,const T lambda2=0, constraint_type mode = PENALTY, const bool pos = false,
122  const int numThreads = -1, Matrix<T>* path = NULL, const int length_path=-1);
123 
124 template <typename T>
125 void lasso2(const Data<T>& X, const AbstractMatrix<T>& G, const AbstractMatrix<T>& DtX,
126  SpMatrix<T>& spalpha,
127  int L, const T constraint, constraint_type mode = PENALTY, const bool pos = false,
128  const int numThreads = -1, Matrix<T>* path = NULL, const int length_path=-1);
129 
131 template <typename T>
132 void lasso_mask(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha, const Matrix<bool>& mask,
133  int L, const T constraint,const T lambda2=0, constraint_type mode = PENALTY, const bool pos = false,
134  const int numThreads = -1);
135 
137 template <typename T>
138 void lassoReweighted(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
139  int L, const T constraint, constraint_type mode, const bool pos,
140  const T sigma,
141  const int numThreads = -1);
142 
144 template <typename T>
145 void coreLARS(Vector<T>& Rdn, Vector<T>& Xdn, Vector<T>& A,
146  Vector<T>& u, Vector<T>& sig,
147  Vector<T>& av, Vector<T>& RUn, Matrix<T>& Un,
148  Matrix<T>& Unds, Matrix<T>& Gs,
149  Matrix<T>& Gsa, Matrix<T>& workT, Matrix<T>& R,
150  const AbstractMatrix<T>& G,T& normX,
151  Vector<int>& ind,Vector<T>& coeffs,const T constraint,
152  const bool ols = false,
153  const bool pos =false,
154  constraint_type mode = L1COEFFS,
155  T* path = NULL, int length_path=-1);
156 
157 template <typename T>
158 void coreLARS2(Vector<T>& DtR, const AbstractMatrix<T>& G,
159  Matrix<T>& Gs,
160  Matrix<T>& Ga,
161  Matrix<T>& invGs,
162  Vector<T>& u,
163  Vector<T>& coeffs,
164  Vector<int>& ind,
165  Matrix<T>& work,
166  T& normX,
167  const constraint_type mode,
168  const T constraint, const bool pos = false,
169  T* pr_path = NULL, int length_path = -1);
170 
171 template <typename T>
172 void coreLARS2W(Vector<T>& DtR, AbstractMatrix<T>& G,
173  Matrix<T>& Gs,
174  Matrix<T>& Ga,
175  Matrix<T>& invGs,
176  Vector<T>& u,
177  Vector<T>& coeffs,
178  const Vector<T>& weights,
179  Vector<int>& ind,
180  Matrix<T>& work,
181  T& normX,
182  const constraint_type mode,
183  const T constraint, const bool pos = false);
184 
186 template <typename T>
187 void downDateLasso(int& j,int& minBasis,T& normX,const bool ols,
188  const bool pos, Vector<T>& Rdn, int* ind,
189  T* coeffs, Vector<T>& sig, Vector<T>& av,
190  Vector<T>& Xdn, Vector<T>& RUn,Matrix<T>& Unm, Matrix<T>& Gsm,
191  Matrix<T>& Gsam, Matrix<T>& Undsm, Matrix<T>& Rm);
192 
193 
194 /* ************************
195  * Iterative thresholding
196  * ************************/
197 
203 template <typename T>
204 void ist(const Matrix<T>& X, const Matrix<T>& D,
205  SpMatrix<T>& spalpha, T lambda, constraint_type mode,
206  const int itermax=500,
207  const T tol = 0.5, const int numThreads = -1);
208 template <typename T>
209 void ist(const Matrix<T>& X, const Matrix<T>& D,
210  Matrix<T>& spalpha, T lambda, constraint_type mode,
211  const int itermax=500,
212  const T tol = 0.5, const int numThreads=-1);
213 
214 
216 template <typename T>
217 void coreIST(const AbstractMatrix<T>& G, Vector<T>& DtR, Vector<T>& coeffs,
218  const T thrs, const int itermax = 500,
219  const T tol = 0.5);
220 
222 template <typename T>
223 void coreISTconstrained(const AbstractMatrix<T>& G, Vector<T>& DtR, Vector<T>& coeffs,
224  const T normX2,
225  const T thrs, const int itermax = 500,
226  const T tol = 0.5);
227 
229 template <typename T>
230 void ist_groupLasso(const Matrix<T>* XT, const Matrix<T>& D,
231  Matrix<T>* alphaT, const int Ngroups,
232  const T lambda, const constraint_type mode,
233  const int itermax = 500,
234  const T tol = 0.5, const int numThreads = -1);
235 
237 template <typename T>
238 void coreGroupIST(const Matrix<T>& G, Matrix<T>& RtD,
239  Matrix<T>& alphat,
240  const T thrs,
241  const int itermax=500,
242  const T tol = 0.5);
243 
244 
246 template <typename T>
247 void coreGroupISTConstrained(const Matrix<T>& G, Matrix<T>& RtD,
248  Matrix<T>& alphat, const T normR,
249  const T eps,
250  const int itermax=500,
251  const T tol = 0.5);
252 
254 template <typename T>
255 T computeError(const T normX2,const Vector<T>& norms,
256  const Matrix<T>& G,const Matrix<T>& RtD,const Matrix<T>& alphat);
257 
259 template <typename T>
260 T computeError(const T normX2,
261  const Matrix<T>& G,const Vector<T>& DtR,const Vector<T>& coeffs,
262  SpVector<T>& coeffs_tmp);
263 
264 /* ******************
265  * Simultaneous OMP
266  * *****************/
267 template <typename T>
268 void somp(const Matrix<T>* X, const Matrix<T>& D, SpMatrix<T>* spalpha,
269  const int Ngroups, const int L, const T* pr_eps, const bool adapt=false,
270  const int numThreads=-1);
271 
272 template <typename T>
273 void somp(const Matrix<T>* X, const Matrix<T>& D, SpMatrix<T>* spalpha,
274  const int Ngroups, const int L, const T eps, const int numThreads=-1);
275 
276 
277 template <typename T>
278 void coreSOMP(const Matrix<T>& X, const Matrix<T>& D, const Matrix<T>& G,
279  Matrix<T>& vM,
280  Vector<int>& rv, const int L, const T eps);
281 
282 /* *********************
283  * Implementation of OMP
284  * *********************/
285 
296 
297 template <typename T>
298 void omp(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
299  const int* pL, const T* peps, const T* pLambda,
300  const bool vecL, const bool vecEps,
301  const bool vecLambda, const int numThreads, Matrix<T>* path) {
302  int L;
303  if (!vecL) {
304  L=*pL;
305  } else {
306  Vector<int> vL(const_cast<int*>(pL),X.n());
307  L=vL.maxval();
308  }
309  spalpha.clear();
310  if (L <= 0) return;
311  const int M = X.n();
312  const int K = D.n();
313  L = MIN(X.m(),MIN(L,K));
314  Matrix<T> vM(L,M);
315  Matrix<int> rM(L,M);
316 
317  ProdMatrix<T> G(D, K < 25000 && M > 10);
318 
319  int NUM_THREADS=init_omp(numThreads);
320 
321  Vector<T>* scoresT=new Vector<T>[NUM_THREADS];
322  Vector<T>* normT=new Vector<T>[NUM_THREADS];
323  Vector<T>* tmpT=new Vector<T>[NUM_THREADS];
324  Vector<T>* RdnT=new Vector<T>[NUM_THREADS];
325  Matrix<T>* UnT=new Matrix<T>[NUM_THREADS];
326  Matrix<T>* UndnT=new Matrix<T>[NUM_THREADS];
327  Matrix<T>* UndsT=new Matrix<T>[NUM_THREADS];
328  Matrix<T>* GsT=new Matrix<T>[NUM_THREADS];
329  for (int i = 0; i<NUM_THREADS; ++i) {
330  scoresT[i].resize(K);
331  normT[i].resize(K);
332  tmpT[i].resize(K);
333  RdnT[i].resize(K);
334  UnT[i].resize(L,L);
335  UnT[i].setZeros();
336  UndnT[i].resize(K,L);
337  UndsT[i].resize(L,L);
338  GsT[i].resize(K,L);
339  }
340 
341  int i;
342 #pragma omp parallel for private(i)
343  for (i = 0; i< M; ++i) {
344 #ifdef _OPENMP
345  int numT=omp_get_thread_num();
346 #else
347  int numT=0;
348 #endif
349  Vector<T> Xi;
350  X.refCol(i,Xi);
351  T normX = Xi.nrm2sq();
352 
353  Vector<int> ind;
354  rM.refCol(i,ind);
355  ind.set(-1);
356 
357  Vector<T> RUn;
358  vM.refCol(i,RUn);
359 
360  Vector<T>& Rdn=RdnT[numT];
361  D.multTrans(Xi,Rdn);
362  coreORMP(scoresT[numT],normT[numT],tmpT[numT],UnT[numT],UndnT[numT],UndsT[numT],
363  GsT[numT],Rdn,G,ind,RUn, normX, vecEps ? peps+i : peps,
364  vecL ? pL+i : pL, vecLambda ? pLambda+i : pLambda,
365  path && i==0 ? path->rawX() : NULL);
366  }
367 
368  delete[](scoresT);
369  delete[](normT);
370  delete[](tmpT);
371  delete[](RdnT);
372  delete[](UnT);
373  delete[](UndnT);
374  delete[](UndsT);
375  delete[](GsT);
376 
378  spalpha.convert(vM,rM,K);
379 };
380 
381 template <typename T>
382 void omp_mask(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha, const Matrix<bool>& mask,
383  const int *pL, const T* peps, const T* pLambda, const bool vecL,
384  const bool vecEps, const bool vecLambda, const int numThreads,
385  Matrix<T>* path) {
386  int L;
387  if (!vecL) {
388  L=*pL;
389  } else {
390  Vector<int> vL(const_cast<int*>(pL),X.n());
391  L=vL.maxval();
392  }
393  spalpha.clear();
394  if (L <= 0) return;
395  const int M = X.n();
396  const int K = D.n();
397  L = MIN(X.m(),MIN(L,K));
398  Matrix<T> vM(L,M);
399  Matrix<int> rM(L,M);
400 
401  ProdMatrix<T> G(D, K < 25000 && M > 10);
402 
403  int NUM_THREADS=init_omp(numThreads);
404 
405  Vector<T>* scoresT=new Vector<T>[NUM_THREADS];
406  Vector<T>* normT=new Vector<T>[NUM_THREADS];
407  Vector<T>* tmpT=new Vector<T>[NUM_THREADS];
408  Vector<T>* RdnT=new Vector<T>[NUM_THREADS];
409  Matrix<T>* UnT=new Matrix<T>[NUM_THREADS];
410  Matrix<T>* UndnT=new Matrix<T>[NUM_THREADS];
411  Matrix<T>* UndsT=new Matrix<T>[NUM_THREADS];
412  Matrix<T>* GsT=new Matrix<T>[NUM_THREADS];
413  ProdMatrix<T>* GT=new ProdMatrix<T>[NUM_THREADS];
414  Matrix<T>* DmaskT=new Matrix<T>[NUM_THREADS];
415  Vector<T>* XmaskT=new Vector<T>[NUM_THREADS];
416  for (int i = 0; i<NUM_THREADS; ++i) {
417  DmaskT[i].resize(D.m(),D.n());
418  XmaskT[i].resize(X.m());
419  scoresT[i].resize(K);
420  normT[i].resize(K);
421  tmpT[i].resize(K);
422  RdnT[i].resize(K);
423  UnT[i].resize(L,L);
424  UnT[i].setZeros();
425  UndnT[i].resize(K,L);
426  UndsT[i].resize(L,L);
427  GsT[i].resize(K,L);
428  }
429 
430  int i;
431 #pragma omp parallel for private(i)
432  for (i = 0; i< M; ++i) {
433 #ifdef _OPENMP
434  int numT=omp_get_thread_num();
435 #else
436  int numT=0;
437 #endif
438  Vector<T> Xi;
439  X.refCol(i,Xi);
440 
441  Vector<int> ind;
442  rM.refCol(i,ind);
443  ind.set(-1);
444 
445  Vector<T> RUn;
446  vM.refCol(i,RUn);
447 
448  Vector<bool> maski;
449  mask.refCol(i,maski);
450  Vector<T>& Rdn=RdnT[numT];
451  if (maski.allfalse()) continue;
452  if (maski.alltrue()) {
453  D.multTrans(Xi,Rdn);
454  T normX = Xi.nrm2sq();
455  coreORMP(scoresT[numT],normT[numT],tmpT[numT],UnT[numT],UndnT[numT],UndsT[numT],
456  GsT[numT],Rdn,G,ind,RUn, normX, vecEps ? peps+i : peps,
457  vecL ? pL+i : pL, vecLambda ? pLambda+i : pLambda,
458  path && i==0 ? path->rawX() : NULL);
459  } else {
460  D.copyMask(DmaskT[numT],maski);
461  Xi.copyMask(XmaskT[numT],maski);
462  T normX = XmaskT[numT].nrm2sq();
463  DmaskT[numT].multTrans(XmaskT[numT],Rdn);
464  GT[numT].setMatrices(DmaskT[numT],false);
465  GT[numT].addDiag(T(1e-10));
466  T eps_mask= (vecEps ? *(peps+i) : *peps)*XmaskT[numT].n()/Xi.n();
467  coreORMP(scoresT[numT],normT[numT],tmpT[numT],
468  UnT[numT],UndnT[numT],UndsT[numT],
469  GsT[numT],Rdn,GT[numT],ind,RUn,
470  normX, &eps_mask, vecL ? pL+i : pL,
471  vecLambda ? pLambda+i : pLambda,
472  path && i==0 ? path->rawX() : NULL);
473 
474  DmaskT[numT].setm(D.m());
475  DmaskT[numT].setn(D.n());
476  XmaskT[numT].setn(X.m());
477  }
478  }
479 
480  delete[](GT);
481  delete[](XmaskT);
482  delete[](DmaskT);
483  delete[](scoresT);
484  delete[](normT);
485  delete[](tmpT);
486  delete[](RdnT);
487  delete[](UnT);
488  delete[](UndnT);
489  delete[](UndsT);
490  delete[](GsT);
491 
493  spalpha.convert(vM,rM,K);
494 };
495 
497 template <typename T>
498 void coreORMPB(Vector<T>& RtD, const AbstractMatrix<T>& G, Vector<int>& ind,
499  Vector<T>& coeffs, T& normX, const int L, const T eps, const T lambda) {
500  const int K = G.n();
501  Vector<T> scores(K);
502  Vector<T> norm(K);
503  Vector<T> tmp(K);
504  Matrix<T> Un(L,L);
505  Matrix<T> Undn(K,L);
506  Matrix<T> Unds(L,L);
507  Matrix<T> Gs(K,L);
508  ind.set(-1);
509  coreORMP(scores,norm,tmp,Un,Undn,Unds,Gs,RtD,G,ind,coeffs,normX,&eps,&L,&lambda);
510 };
511 
513 template <typename T>
514 void coreORMP(Vector<T>& scores, Vector<T>& norm, Vector<T>& tmp, Matrix<T>& Un,
515  Matrix<T>& Undn, Matrix<T>& Unds, Matrix<T>& Gs, Vector<T>& Rdn,
516  const AbstractMatrix<T>& G,
517  Vector<int>& ind, Vector<T>& RUn,
518  T& normX, const T* peps, const int* pL, const T* plambda,
519  T* path) {
520  const T eps = abs<T>(*peps);
521  const int L = MIN(*pL,Gs.n());
522  const T lambda=*plambda;
523  if ((normX <= eps) || L == 0) return;
524  const int K = scores.n();
525  scores.copy(Rdn);
526  norm.set(T(1.0));
527  Un.setZeros();
528 
529  // permit unsafe low level access
530  T* const prUn = Un.rawX();
531  T* const prUnds = Unds.rawX();
532  T* const prUndn = Undn.rawX();
533  T* const prGs = Gs.rawX();
534  T* const prRUn= RUn.rawX();
535  if (path)
536  memset(path,0,K*L*sizeof(T));
537 
538  int j;
539  for (j = 0; j<L; ++j) {
540  const int currentInd=scores.fmax();
541  if (norm[currentInd] < 1e-8) {
542  ind[j]=-1;
543  break;
544  }
545  const T invNorm=T(1.0)/sqrt(norm[currentInd]);
546  const T RU=Rdn[currentInd]*invNorm;
547  const T delta = RU*RU;
548  if (delta < 2*lambda) {
549  break;
550  }
551 
552  RUn[j]=RU;
553  normX -= delta;
554  ind[j]=currentInd;
555  //for (int k = 0; k<j; ++k) prUn[j*L+k]=0.0;
556  //prUn[j*L+j]=T(1.0);
557 
558  // for (int k = 0; k<j; ++k) prUnds[k*L+j]=prUndn[k*K+currentInd];
559  // MGS algorithm, Update Un
560  // int iter = norm[currentInd] < 0.5 ? 2 : 1;
561  //int iter=1;
562  // for (int k = 0; k<iter; ++k) {
564  // T scal=-cblas_dot<T>(j+1-l,prUn+j*L+l,1,prUnds+l*L+l,1);
565  // T scal = -prUnds[l*L+j];
566  // cblas_axpy<T>(l+1,scal,prUn+l*L,1,prUn+j*L,1);
567  // }
568  // }
569 
570  prUn[j*L+j]=-T(1.0);
571  cblas_copy<T>(j,prUndn+currentInd,K,prUn+j*L,1);
572  cblas_trmv<T>(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit,j,prUn,L,prUn+j*L,1);
573  cblas_scal<T>(j+1,-invNorm,prUn+j*L,1);
574 
575  if (j == L-1 || (normX <= eps)) {
576  ++j;
577  break;
578  }
579 
580  if (path) {
581  T* last_path=path+(L-1)*K;
582  cblas_copy<T>(j+1,prRUn,1,last_path,1);
583  cblas_trmv<T>(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit,
584  j+1,prUn,L,last_path,1);
585  for (int k = 0; k<=j; ++k) {
586  path[j*K+ind[k]]=last_path[k];
587  }
588  }
589 
590  // update the variables Gs, Undn, Unds, Rdn, norm, scores
591  Vector<T> Gsj;
592  Gs.refCol(j,Gsj);
593  G.copyCol(currentInd,Gsj);
594  cblas_gemv<T>(CblasColMajor,CblasNoTrans,K,j+1,T(1.0),prGs,K,prUn+j*L,1,
595  T(0.0),prUndn+j*K,1);
596  // prUnds[j*L+j] = prUndn[j*K+currentInd];
597  Vector<T> Undnj;
598  Undn.refCol(j,Undnj);
599  Rdn.add(Undnj,-RUn[j]);
600  tmp.sqr(Undnj);
601  norm.sub(tmp);
602  scores.sqr(Rdn);
603  scores.div(norm);
604  for (int k = 0; k<=j; ++k) scores[ind[k]]=T();
605  }
606  // compute the final coefficients
607  cblas_trmv<T>(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit,
608  j,prUn,L,prRUn,1);
609  if (path) {
610  memset(path+(L-1)*K,0,L*sizeof(T));
611  for (int k = 0; k<j; ++k) {
612  path[(j-1)*K+ind[k]]=prRUn[k];
613  }
614  }
615 };
616 
617 
618 
619 /* **************
620  * LARS - Lasso
621  * **************/
622 
637 
638 template <typename T>
639 void lasso(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
640  int L, const T lambda, const T lambda2, constraint_type mode,
641  const bool pos, const bool ols, const int numThreads,
642  Matrix<T>* path, const int length_path) {
643  ProdMatrix<T> G(D, X.n() > 10 && D.n() < 50000);
644  G.addDiag(MAX(lambda2,1e-10));
645  ProdMatrix<T> DtX(D,X,false);
646  lasso(X,G,DtX,spalpha,L,lambda,mode,pos,ols,numThreads,path,length_path);
647 }
648 
649 template <typename T>
650 void lasso(const Data<T>& X, const AbstractMatrix<T>& G,
651  const AbstractMatrix<T>& DtX, SpMatrix<T>& spalpha,
652  int L, const T lambda, constraint_type mode,
653  const bool pos, const bool ols, const int numThreads,
654  Matrix<T>* path, const int length_path) {
655 
656  spalpha.clear();
657  const int M = X.n();
658  const int K = G.n();
659  Matrix<T> vM;
660  Matrix<int> rM;
661  vM.resize(L,M);
662  rM.resize(L,M);
663 
664  if (L <= 0) return;
665  if (path) path->setZeros();
666 
667  int NUM_THREADS=init_omp(numThreads);
668 
669  //ProdMatrix<T> G(D, K < 25000 && M > 10);
670 
671  Vector<T>* RdnT=new Vector<T>[NUM_THREADS];
672  Vector<T>* XdnT =new Vector<T>[NUM_THREADS];
673  Vector<T>* AT=new Vector<T>[NUM_THREADS];
674  Vector<T>* uT=new Vector<T>[NUM_THREADS];
675  Vector<T>* sigT=new Vector<T>[NUM_THREADS];
676  Vector<T>* avT=new Vector<T>[NUM_THREADS];
677  Vector<T>* RUnT = new Vector<T>[NUM_THREADS];
678  Matrix<T>* UnT=new Matrix<T>[NUM_THREADS];
679  Matrix<T>* RT=new Matrix<T>[NUM_THREADS];
680  Matrix<T>* UndsT=new Matrix<T>[NUM_THREADS];
681  Matrix<T>* GsT=new Matrix<T>[NUM_THREADS];
682  Matrix<T>* GsaT=new Matrix<T>[NUM_THREADS];
683  Matrix<T>* workT=new Matrix<T>[NUM_THREADS];
684  for (int i = 0; i<NUM_THREADS; ++i) {
685  RdnT[i].resize(K);
686  if (ols) XdnT[i].resize(K);
687  AT[i].resize(K);
688  uT[i].resize(L);
689  sigT[i].resize(L);
690  avT[i].resize(L);
691  if (ols) RUnT[i].resize(L);
692  UnT[i].resize(L,L);
693  UnT[i].setZeros();
694  UndsT[i].resize(L,L);
695  UndsT[i].setZeros();
696  GsT[i].resize(K,L);
697  GsaT[i].resize(L,L);
698  workT[i].resize(K,2);
699  RT[i].resize(L,L);
700  }
701 
702  Vector<T> norms;
703  X.norm_2sq_cols(norms);
704  int i;
705 #pragma omp parallel for private(i)
706  for (i = 0; i< M; ++i) {
707 #ifdef _OPENMP
708  int numT=omp_get_thread_num();
709 #else
710  int numT=0;
711 #endif
712  T normX = norms[i];
713 
714  Vector<int> ind;
715  rM.refCol(i,ind);
716  Vector<T> coeffs;
717  vM.refCol(i,coeffs);
718  coeffs.setZeros();
719 
720  Vector<T>& Rdn=RdnT[numT];
721  DtX.copyCol(i,Rdn);
722  coreLARS(Rdn,XdnT[numT], AT[numT], uT[numT], sigT[numT], avT[numT],
723  RUnT[numT], UnT[numT], UndsT[numT], GsT[numT], GsaT[numT],
724  workT[numT],RT[numT],G,normX, ind,coeffs,lambda,ols,pos,
725  mode,path && i==0 ? path->rawX() : NULL, length_path);
726  }
727 
728  delete[](RdnT);
729  delete[](XdnT);
730  delete[](AT);
731  delete[](uT);
732  delete[](sigT);
733  delete[](avT);
734  delete[](RUnT);
735  delete[](UnT);
736  delete[](RT);
737  delete[](UndsT);
738  delete[](GsT);
739  delete[](GsaT);
740  delete[](workT);
741 
743  spalpha.convert(vM,rM,K);
744 };
745 
747 template <typename T>
748 void coreLARS(Vector<T>& Rdnv, Vector<T>& Xdnv, Vector<T>& Av,
749  Vector<T>& uv, Vector<T>& sigv, Vector<T>& avv, Vector<T>& RUnv,
750  Matrix<T>& Unm, Matrix<T>& Undsm, Matrix<T>& Gsm,
751  Matrix<T>& Gsam, Matrix<T>& workm, Matrix<T>& Rm,
752  const AbstractMatrix<T>& Gm,T& normX,
753  Vector<int>& indv,Vector<T>& coeffsv,const T constraint,
754  const bool ols,const bool pos, constraint_type mode,
755  T* path, int length_path) {
756  if (mode == L2ERROR && normX < constraint) return;
757 
758  const int LL = Gsm.n();
759  const int K = Gsm.m();
760  const int L = MIN(LL,K);
761  if (length_path <= 1) length_path=4*L;
762  // permit unsafe fast low level access
763  T* const Rdn = Rdnv.rawX();
764  T* const Xdn = Xdnv.rawX();
765  T* const A = Av.rawX();
766  T* const u = uv.rawX();
767  T* const sig = sigv.rawX();
768  T* const av = avv.rawX();
769  T* const RUn = RUnv.rawX();
770  T* const Un = Unm.rawX();
771  T* const Unds = Undsm.rawX();
772  T* const Gs = Gsm.rawX();
773  T* const Gsa = Gsam.rawX();
774  T* const work = workm.rawX();
775  //T* const G = Gm.rawX();
776  T* const R = Rm.rawX();
777  int* ind = indv.rawX();
778  T* coeffs = coeffsv.rawX();
779 
780  coeffsv.setZeros();
781  indv.set(-1);
782 
783  if (ols) Xdnv.copy(Rdnv);
784  int currentInd= pos ? Rdnv.max() : Rdnv.fmax();
785  bool newAtom=true;
786  T Cmax;
787  int iter=1;
788  T thrs = 0.0;
789 
790  int* const ind_orig = ind;
791  T* const coeffs_orig = coeffs;
792 
793  int j;
794  for (j = 0; j<L; ++j) {
795  if (newAtom) {
796  ind[j]=currentInd;
797 
798  if (pos) {
799  Cmax = Rdn[currentInd];
800  sig[j]=1.0;
801  } else {
802  Cmax = abs<T>(Rdn[currentInd]);
803  sig[j] = SIGN(Rdn[currentInd]);
804  }
805  for (int k = 0; k<=j; ++k) Un[j*L+k]=0.0;
806  Un[j*L+j]=1.0;
807  Gm.extract_rawCol(currentInd,Gs+K*j);
808  for (int k = 0; k<j; ++k) Gs[K*j+ind[k]] *= sig[k];
809  if (sig[j] < 0) {
810  Rdn[currentInd]=-Rdn[currentInd];
811  if (ols) Xdn[currentInd]=-Xdn[currentInd];
812  cblas_scal<T>(K,sig[j],Gs+K*j,1);
813  cblas_scal<T>(j+1,sig[j],Gs+currentInd,K);
814  }
815  cblas_copy<T>(j+1,Gs+currentInd,K,Gsa+j*L,1);
816  for (int k = 0; k<j; ++k) Gsa[k*L+j]=Gsa[j*L+k];
817 
818  // <d_j,d_i>
819  cblas_copy<T>(j,Gsa+j*L,1,Unds+j,L);
820  // <U_j final,d_i>
821  cblas_trmv<T>(CblasColMajor,CblasUpper,CblasTrans,CblasNonUnit,
822  j+1,Un,L,Unds+j,L);
823  // norm2
824  T norm2=Gsa[j*L+j];
825  for (int k = 0; k<j; ++k) norm2 -= Unds[k*L+j]*Unds[k*L+j];
826  if (norm2 < 1e-15) {
827  ind[j]=-1;
828  // cerr << "bad exit" << endl;
829  break;
830  }
831 
832  // int iter2 = norm2 < 0.5 ? 2 : 1;
833  // for(int k = 0; k<iter2; ++k) {
834  // for (int l = 0; l<j; ++l) {
835  // T scal=-cblas_dot<T>(j+1-l,Un+j*L+l,1,Unds+l*L+l,1);
836  // cblas_axpy<T>(l+1,scal,Un+l*L,1,Un+j*L,1);
837  // }
838  // }
839  Un[j*L+j]=-T(1.0);
840  cblas_copy<T>(j,Unds+j,L,Un+j*L,1);
841  cblas_trmv<T>(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit,j,Un,L,Un+j*L,1);
842 
844  T invNorm=1.0/sqrt(norm2);
845  cblas_scal<T>(j+1,-invNorm,Un+j*L,1);
846  Unds[j*L+j]=cblas_dot<T>(j+1,Un+j*L,1,Gsa+j*L,1);
847  }
848 
849  for (int k = 0; k<=j; ++k) u[k]=T(1.0);
850  cblas_trmv<T>(CblasColMajor,CblasUpper,CblasTrans,CblasNonUnit,
851  j+1,Un,L,u,1);
852 
853  T a = T(1.0)/cblas_nrm2<T>(j+1,u,1);
854 
855  cblas_trmv<T>(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit,
856  j+1,Un,L,u,1);
857  cblas_scal<T>(j+1,a,u,1);
858 
859  cblas_gemv<T>(CblasColMajor,CblasNoTrans,K,j+1,T(1.0),Gs,K,u,1,T(0.0),A,1);
860 
861  T potentNorm=0.0;
862  if (!ols) {
863  for (int k = 0; k<=j; ++k) potentNorm += Rdn[ind[k]]*u[k];
864  }
865 
866  if (pos) {
867  for (int k = 0; k<K; ++k) {
868  T diff = a-A[k];
869  work[k]= diff <= 0 ? INFINITY : (Cmax-Rdn[k])/diff;
870  }
871  for (int k = 0; k<=j; ++k) {
872  work[ind[k]]=INFINITY;
873  }
874  for (int k = 0; k<K; ++k)
875  if (work[k] <=0) work[k]=INFINITY;
876  currentInd =cblas_iamin<T>(K,work,1);
877  } else {
878  memset(work,0,2*K*sizeof(T));
879  for (int k = 0; k<=j; ++k) {
880  const int index=2*ind[k];
881  work[index]=INFINITY;
882  work[index+1]=INFINITY;
883  }
884  for (int k = 0; k<K; ++k) {
885  const int index=2*k;
886  if (!work[index]) {
887  const T diff1=a-A[k];
888  work[index]= diff1 <= 0 ? INFINITY : (Cmax-Rdn[k])/diff1;
889  const T diff2=a+A[k];
890  work[index+1]=diff2 <= 0 ? INFINITY : (Cmax+Rdn[k])/diff2;
891  }
892  }
893  currentInd =cblas_iamin<T>(2*K,work,1);
894  }
895  T gamma=work[currentInd];
896  T gammaMin=0;
897  int minBasis=0;
898 
899  //if (j == L-1) gamma=potentNorm;
900 
901  if (mode == PENALTY) {
902  gamma=MIN(gamma,(Cmax-constraint)/a);
903  }
904 
905 // if (j > 0) {
906  vDiv<T>(j+1,coeffs,u,work);
907  cblas_scal<T>(j+1,-T(1.0),work,1);
909  for (int k=0; k<=j; ++k)
910  if (coeffs[k]==0 || work[k] <=0) work[k]=INFINITY;
911  minBasis=cblas_iamin<T>(j+1,work,1);
912  gammaMin=work[minBasis];
913  if (gammaMin < gamma) gamma=gammaMin;
914  // }
915 
916  if (mode == L1COEFFS) {
917  T Tu = 0.0;
918  for (int k = 0; k<=j; ++k) Tu += u[k];
919 
920  if (Tu > EPSILON)
921  gamma= MIN(gamma,(constraint-thrs)/Tu);
922  thrs+=gamma*Tu;
923  }
924 
925  // compute the norm of the residdual
926 
927  if (ols == 0) {
928  const T t = gamma*gamma - 2*gamma*potentNorm;
929  if (t > 0 || std::isnan(t) || std::isinf(t)) {
930  // cerr << "bad bad exit" << endl;
931  // cerr << t << endl;
932  ind[j]=-1;
933  break;
934  }
935  normX += t;
936  } else {
937  // plan the last orthogonal projection
938  if (newAtom) {
939  RUn[j]=0.0;
940  for (int k = 0; k<=j; ++k) RUn[j] += Xdn[ind[k]]*
941  Un[j*L+k];
942  normX -= RUn[j]*RUn[j];
943  }
944  }
945 
946  // Update the coefficients
947  cblas_axpy<T>(j+1,gamma,u,1,coeffs,1);
948 
949  if (pos) {
950  for (int k = 0; k<j+1; ++k)
951  if (coeffs[k] < 0) coeffs[k]=0;
952  }
953 
954  cblas_axpy<T>(K,-gamma,A,1,Rdn,1);
955  if (!pos) currentInd/= 2;
956  if (path) {
957  for (int k = 0; k<=j; ++k)
958  path[iter*K+ind[k]]=coeffs[k]*sig[k];
959  }
960 
961  if (gamma == gammaMin) {
962  downDateLasso<T>(j,minBasis,normX,ols,pos,Rdnv,ind,coeffs,sigv,
963  avv,Xdnv, RUnv, Unm, Gsm, Gsam,Undsm,Rm);
964  newAtom=false;
965  Cmax=abs<T>(Rdn[ind[0]]);
966  --j;
967  } else {
968  newAtom=true;
969  }
970  ++iter;
971 
972  if (mode == PENALTY) {
973  thrs=abs<T>(Rdn[ind[0]]);
974  }
975 
976  if ((j == L-1) ||
977  (mode == PENALTY && (thrs - constraint < 1e-15)) ||
978  (mode == L1COEFFS && (thrs - constraint > -1e-15)) ||
979  (newAtom && mode == L2ERROR && (normX - constraint < 1e-15)) ||
980  (normX < 1e-15) ||
981  (iter >= length_path)) {
982  // cerr << "exit" << endl;
983  // PRINT_F(thrs)
984  // PRINT_F(constraint)
985  // PRINT_F(normX)
986  break;
987  }
988 
989  }
990  if (ols) {
991  cblas_copy<T>(j+1,RUn,1,coeffs,1);
992  cblas_trmv<T>(CblasColMajor,CblasUpper,CblasNoTrans,CblasNonUnit,
993  j+1,Un,L,coeffs,1);
994  }
995  vMul<T>(j+1,coeffs,sig,coeffs);
996 };
997 
999 template <typename T>
1000 inline void downDateLasso(int& j,int& minBasis,T& normX,const bool ols,
1001  const bool pos,
1002  Vector<T>& Rdnv, int* ind,
1003  T* coeffs, Vector<T>& sigv, Vector<T>& avv,
1004  Vector<T>& Xdnv, Vector<T>& RUnv,Matrix<T>& Unm, Matrix<T>& Gsm,
1005  Matrix<T>& Gsam, Matrix<T>& Undsm, Matrix<T>& Rm) {
1006  int k,l;
1007  const int L = Gsm.n();
1008  const int K = Gsm.m();
1009  T* const Rdn = Rdnv.rawX();
1010  T* const Xdn = Xdnv.rawX();
1011  T* const sig = sigv.rawX();
1012  T* const av = avv.rawX();
1013  T* const RUn = RUnv.rawX();
1014  T* const Un = Unm.rawX();
1015  T* const Unds = Undsm.rawX();
1016  T* const Gs = Gsm.rawX();
1017  T* const Gsa = Gsam.rawX();
1018  T* const R = Rm.rawX();
1019 
1020  int indB=ind[minBasis];
1021 
1022  if (!pos && sig[minBasis] < 0) {
1023  // Update Rdn
1024  Rdn[indB]=-Rdn[indB];
1025  if (ols) Xdn[indB]=-Xdn[indB];
1026  }
1027 
1028  int num=j-minBasis;
1029  for (int k = 0; k<num*num;++k) R[k]=0.0;
1030  for (int k = 0; k<num; ++k) R[k*num+k]=1.0;
1031  // Update Un
1032  for (int k = minBasis+1; k<=j; ++k) {
1033  T a = -Un[k*L+minBasis]/Un[minBasis*L+minBasis];
1034  av[k-minBasis-1] = a;
1035  cblas_axpy<T>(minBasis,a,Un+minBasis*L,1,Un+k*L,1);
1036  }
1037  for (int k = minBasis+1; k<=j; ++k) {
1038  cblas_copy<T>(minBasis,Un+k*L,1,Un+(k-1)*L,1);
1039  cblas_copy<T>(num,Un+k*L+minBasis+1,1,Un+(k-1)*L+minBasis,1);
1040  }
1041  T alpha=1.0;
1042  T alphab,gamma,lambda;
1043  for (int k = 0; k<num; ++k) {
1044  alphab=alpha+av[k]*av[k];
1045  R[k*num+k]=sqrt(alphab/alpha);
1046  gamma=av[k]*R[k*num+k]/alphab;
1047  alpha=alphab;
1048  cblas_copy<T>(num-k-1,av+k+1,1,R+k*num+k+1,1);
1049  cblas_scal<T>(num-k-1,gamma,R+k*num+k+1,1);
1050  }
1051  if (num > 0) {
1052  trtri<T>(low,nonUnit,num,R,num);
1053  cblas_trmm<T>(CblasColMajor,CblasRight,CblasLower,CblasTrans,CblasNonUnit,
1054  j,num,T(1.0),R,num,Un+minBasis*L,L);
1055  }
1056 
1057  // Update Unds
1058  for (int k = minBasis+1; k<=j; ++k)
1059  cblas_axpy<T>(j-minBasis,av[k-minBasis-1],Unds+minBasis*L+minBasis+1,1,
1060  Unds+k*L+minBasis+1,1);
1061  for (int k = 0; k<minBasis; ++k)
1062  for (int l = minBasis+1; l<=j; ++l)
1063  Unds[k*L+l-1]=Unds[k*L+l];
1064  for (int k = minBasis+1; k<=j; ++k)
1065  cblas_copy<T>(j-minBasis,Unds+k*L+minBasis+1,1,Unds+(k-1)*L+minBasis,1);
1066  if (num > 0)
1067  cblas_trmm<T>(CblasColMajor,CblasRight,CblasLower,CblasTrans,CblasNonUnit,
1068  j-minBasis,num,T(1.0),R,num,Unds+minBasis*L+minBasis,L);
1069  for (int k = minBasis+1; k<=j; ++k)
1070  for (int l = 0; l<k; ++l) Unds[k*L+l]=0.0;
1071 
1072  // Update Gs
1073  for (int k = minBasis+1; k<=j; ++k) {
1074  cblas_copy<T>(K,Gs+k*K,1,Gs+(k-1)*K,1);
1075  }
1076  if (!pos && sig[minBasis] < T(0.0)) cblas_scal<T>(j,T(-1.0),Gs+indB,K);
1077  // Update Gsa
1078  for (int k = minBasis+1; k<=j; ++k) {
1079  cblas_copy<T>(minBasis,Gsa+k*L,1,Gsa+(k-1)*L,1);
1080  cblas_copy<T>(j-minBasis,Gsa+k*L+minBasis+1,1,Gsa+(k-1)*L+minBasis,1);
1081  }
1082  for (int k = 0; k<minBasis; ++k) {
1083  for (int l = minBasis+1; l<=j; ++l) Gsa[k*L+l-1]=Gsa[k*L+l];
1084  }
1085 
1086  // Update sig
1087  for (int k = minBasis+1; k<=j && !pos; ++k) sig[k-1]=sig[k];
1088  // Update ind
1089  for (int k = minBasis+1; k<=j; ++k) ind[k-1]=ind[k];
1090  ind[j]=-1;
1091 
1092  for (int k = minBasis+1; k<=j; ++k) coeffs[k-1]=coeffs[k];
1093  coeffs[j]=0.0;
1094 
1095  if (ols) {
1096  // Update RUn and normX
1097  for (int k = minBasis; k<=j; ++k)
1098  normX += RUn[k]*RUn[k];
1099  for (int k = minBasis; k<j; ++k) {
1100  RUn[k]=0.0;
1101  for (int l = 0; l<=k; ++l) RUn[k] += Xdn[ind[l]]*
1102  Un[k*L+l];
1103  normX -= RUn[k]*RUn[k];
1104  }
1105  }
1106 
1107  // Update j
1108  --j;
1109 }
1110 
1112 template <typename T>
1113 void lassoReweighted(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
1114  int L, const T constraint, constraint_type mode, const bool pos,
1115  const T sigma,
1116  const int numThreads) {
1117  spalpha.clear();
1118  const int M = X.n();
1119  const int K = D.n();
1120  Matrix<T> vM;
1121  Matrix<int> rM;
1122  vM.resize(L,M);
1123  rM.resize(L,M);
1124  const int iterR = 30;
1125 
1126  if (L <= 0) return;
1127 
1128  int NUM_THREADS=init_omp(numThreads);
1129 
1130  //ProdMatrix<T> G(D, K < 25000 && M > 10);
1131  ProdMatrix<T> G(D, K < 50000);
1132  //Matrix<T> G;
1133  //D.XtX(G);
1134  G.addDiag(1e-10);
1135 
1136  Vector<T>* DtRT=new Vector<T>[NUM_THREADS];
1137  Vector<T>* DtRRT=new Vector<T>[NUM_THREADS];
1138  Vector<T>* uT=new Vector<T>[NUM_THREADS];
1139  Vector<T>* weightsT=new Vector<T>[NUM_THREADS];
1140  Vector<int>* inddT=new Vector<int>[NUM_THREADS];
1141  Matrix<T>* GsT=new Matrix<T>[NUM_THREADS];
1142  Matrix<T>* GaT=new Matrix<T>[NUM_THREADS];
1143  Matrix<T>* invGsT=new Matrix<T>[NUM_THREADS];
1144  Matrix<T>* workT=new Matrix<T>[NUM_THREADS];
1145  Matrix<T>* GT=new Matrix<T>[NUM_THREADS];
1146  for (int i = 0; i<NUM_THREADS; ++i) {
1147  DtRT[i].resize(K);
1148  DtRRT[i].resize(K);
1149  uT[i].resize(K);
1150  weightsT[i].resize(K);
1151  GT[i].resize(K,K);
1152  inddT[i].resize(K);
1153  GsT[i].resize(L,L);
1154  invGsT[i].resize(L,L);
1155  GaT[i].resize(K,L);
1156  workT[i].resize(K,3);
1157  workT[i].setZeros();
1158  }
1159 
1160  int i;
1161 #pragma omp parallel for private(i)
1162  for (i = 0; i< M; ++i) {
1163 #ifdef _OPENMP
1164  int numT=omp_get_thread_num();
1165 #else
1166  int numT=0;
1167 #endif
1168  Vector<T> Xi;
1169  X.refCol(i,Xi);
1170  T normXo = Xi.nrm2sq();
1171  T normX = normXo;
1172 
1173  Vector<int> ind;
1174  rM.refCol(i,ind);
1175  Vector<T> coeffs;
1176  vM.refCol(i,coeffs);
1177  Vector<T>& DtR=DtRT[numT];
1178  Vector<T>& DtRR = DtRRT[numT];
1179  D.multTrans(Xi,DtR);
1180  DtRR.copy(DtR);
1181  coreLARS2(DtRR,G,GsT[numT],GaT[numT],invGsT[numT],uT[numT],coeffs,
1182  ind,workT[numT],normX,mode,constraint,pos);
1183  //Matrix<T>& GG = GT[numT];
1184  Vector<T>& weights = weightsT[numT];
1185  //Vector<int>& indd = inddT[numT];
1186  for (int j = 0; j<iterR; ++j) {
1187  const T sig = sigma*pow(0.7,iterR-1-j);
1188  weights.set(sig);
1189  for (int k = 0; k<K; ++k) {
1190  if (ind[k] != -1) {
1191  weights[ind[k]] = MAX(1e-4,sig*exp(-sig*abs<T>(coeffs[k])));
1192  } else {
1193  break;
1194  }
1195  }
1196  DtRR.copy(DtR);
1197  normX=normXo;
1198  coreLARS2W(DtRR,G,GsT[numT],GaT[numT],invGsT[numT],uT[numT],coeffs,weights,
1199  ind,workT[numT],normX,mode,constraint,pos);
1200  }
1201  }
1202 
1203  delete[](DtRT);
1204  delete[](DtRRT);
1205  delete[](inddT);
1206  delete[](uT);
1207  delete[](weightsT);
1208  delete[](GsT);
1209  delete[](GT);
1210  delete[](GaT);
1211  delete[](invGsT);
1212  delete[](workT);
1213 
1215  spalpha.convert(vM,rM,K);
1216 
1217 }
1218 
1219 template <typename T>
1220 void lassoWeight(const Matrix<T>& X, const Matrix<T>& D, const Matrix<T>& weights,
1221  SpMatrix<T>& spalpha,
1222  int L, const T constraint, constraint_type mode, const bool pos,
1223  const int numThreads) {
1224 
1225  spalpha.clear();
1226  const int M = X.n();
1227  const int K = D.n();
1228  Matrix<T> vM;
1229  Matrix<int> rM;
1230  vM.resize(L,M);
1231  rM.resize(L,M);
1232 
1233  if (L <= 0) return;
1234 
1235  int NUM_THREADS=init_omp(numThreads);
1236 
1237  //ProdMatrix<T> G(D, K < 25000 && M > 10);
1238  ProdMatrix<T> G(D, K < 50000);
1239  //Matrix<T> G;
1240  //D.XtX(G);
1241  G.addDiag(1e-10);
1242 
1243  Vector<T>* DtRT=new Vector<T>[NUM_THREADS];
1244  Vector<T>* uT=new Vector<T>[NUM_THREADS];
1245  Matrix<T>* GsT=new Matrix<T>[NUM_THREADS];
1246  Matrix<T>* GaT=new Matrix<T>[NUM_THREADS];
1247  Matrix<T>* invGsT=new Matrix<T>[NUM_THREADS];
1248  Matrix<T>* workT=new Matrix<T>[NUM_THREADS];
1249  for (int i = 0; i<NUM_THREADS; ++i) {
1250  DtRT[i].resize(K);
1251  uT[i].resize(K);
1252  uT[i].setZeros();
1253  GsT[i].resize(L,L);
1254  invGsT[i].resize(L,L);
1255  GaT[i].resize(K,L);
1256  workT[i].resize(K,3);
1257  workT[i].setZeros();
1258  }
1259 
1260  int i;
1261 #pragma omp parallel for private(i)
1262  for (i = 0; i< M; ++i) {
1263 #ifdef _OPENMP
1264  int numT=omp_get_thread_num();
1265 #else
1266  int numT=0;
1267 #endif
1268  Vector<T> Xi;
1269  X.refCol(i,Xi);
1270  T normX = Xi.nrm2sq();
1271 
1272  Vector<int> ind;
1273  rM.refCol(i,ind);
1274  Vector<T> coeffs;
1275  vM.refCol(i,coeffs);
1276 
1277  Vector<T>& DtR=DtRT[numT];
1278  D.multTrans(Xi,DtR);
1279  Vector<T> we;
1280  weights.refCol(i,we);
1281 
1282  coreLARS2W(DtR,G,GsT[numT],GaT[numT],invGsT[numT],uT[numT],coeffs,we,
1283  ind,workT[numT],normX,mode,constraint,pos);
1284  }
1285 
1286  delete[](DtRT);
1287  delete[](uT);
1288  delete[](GsT);
1289  delete[](GaT);
1290  delete[](invGsT);
1291  delete[](workT);
1292 
1294  spalpha.convert(vM,rM,K);
1295 };
1296 
1297 template <typename T>
1298 void lassoWeightPreComputed(const Matrix<T>& X, const Matrix<T>& G, const Matrix<T>& DtR, const Matrix<T>& weights,
1299  SpMatrix<T>& spalpha,
1300  int L, const T constraint, constraint_type mode, const bool pos,
1301  const int numThreads) {
1302 
1303  spalpha.clear();
1304  const int M = X.n();
1305  const int K = G.n();
1306  Matrix<T> vM;
1307  Matrix<int> rM;
1308  vM.resize(L,M);
1309  rM.resize(L,M);
1310 
1311  if (L <= 0) return;
1312 
1313  int NUM_THREADS=init_omp(numThreads);
1314 
1315  Vector<T>* DtRT=new Vector<T>[NUM_THREADS];
1316  Vector<T>* uT=new Vector<T>[NUM_THREADS];
1317  Matrix<T>* GsT=new Matrix<T>[NUM_THREADS];
1318  Matrix<T>* GaT=new Matrix<T>[NUM_THREADS];
1319  Matrix<T>* invGsT=new Matrix<T>[NUM_THREADS];
1320  Matrix<T>* workT=new Matrix<T>[NUM_THREADS];
1321  for (int i = 0; i<NUM_THREADS; ++i) {
1322  DtRT[i].resize(K);
1323  uT[i].resize(K);
1324  uT[i].setZeros();
1325  GsT[i].resize(L,L);
1326  invGsT[i].resize(L,L);
1327  GaT[i].resize(K,L);
1328  workT[i].resize(K,3);
1329  workT[i].setZeros();
1330  }
1331 
1332  int i;
1333 #pragma omp parallel for private(i)
1334  for (i = 0; i< M; ++i) {
1335 #ifdef _OPENMP
1336  int numT=omp_get_thread_num();
1337 #else
1338  int numT=0;
1339 #endif
1340  Vector<T> Xi;
1341  X.refCol(i,Xi);
1342  T normX = Xi.nrm2sq();
1343 
1344  Vector<int> ind;
1345  rM.refCol(i,ind);
1346  Vector<T> coeffs;
1347  vM.refCol(i,coeffs);
1348 
1349  Vector<T>& DtRi=DtRT[numT];
1350  DtR.copyCol(i,DtRi);
1351  Vector<T> we;
1352  weights.refCol(i,we);
1353 
1354  coreLARS2W(DtRi,G,GsT[numT],GaT[numT],invGsT[numT],uT[numT],coeffs,we,
1355  ind,workT[numT],normX,mode,constraint,pos);
1356  }
1357 
1358  delete[](DtRT);
1359  delete[](uT);
1360  delete[](GsT);
1361  delete[](GaT);
1362  delete[](invGsT);
1363  delete[](workT);
1364 
1366  spalpha.convert(vM,rM,K);
1367 };
1368 
1370 template <typename T>
1371 void lasso_mask(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha, const Matrix<bool>& mask,
1372  int L, const T constraint,const T lambda2, constraint_type mode, const bool pos,
1373  const int numThreads) {
1374  spalpha.clear();
1375  const int M = X.n();
1376  const int K = D.n();
1377  Matrix<T> vM;
1378  Matrix<int> rM;
1379  vM.resize(L,M);
1380  rM.resize(L,M);
1381 
1382  if (L <= 0) return;
1383 
1384  int NUM_THREADS=init_omp(numThreads);
1385 
1386  ProdMatrix<T> G(D,K < 25000 && M > 10);
1387  G.addDiag(MAX(lambda2,1e-10));
1388 
1389  Vector<T>* DtRT=new Vector<T>[NUM_THREADS];
1390  Vector<T>* uT=new Vector<T>[NUM_THREADS];
1391  Vector<T>* XmaskT=new Vector<T>[NUM_THREADS];
1392  Matrix<T>* GsT=new Matrix<T>[NUM_THREADS];
1393  ProdMatrix<T>* GT=new ProdMatrix<T>[NUM_THREADS];
1394  Matrix<T>* DmaskT=new Matrix<T>[NUM_THREADS];
1395  Matrix<T>* GaT=new Matrix<T>[NUM_THREADS];
1396  Matrix<T>* invGsT=new Matrix<T>[NUM_THREADS];
1397  Matrix<T>* workT=new Matrix<T>[NUM_THREADS];
1398  for (int i = 0; i<NUM_THREADS; ++i) {
1399  DmaskT[i].resize(D.m(),D.n());
1400  DtRT[i].resize(K);
1401  uT[i].resize(K);
1402  XmaskT[i].resize(X.m());
1403  uT[i].setZeros();
1404  GsT[i].resize(L,L);
1405  invGsT[i].resize(L,L);
1406  GaT[i].resize(K,L);
1407  workT[i].resize(K,3);
1408  workT[i].setZeros();
1409  }
1410 
1411  int i;
1412 #pragma omp parallel for private(i)
1413  for (i = 0; i< M; ++i) {
1414 #ifdef _OPENMP
1415  int numT=omp_get_thread_num();
1416 #else
1417  int numT=0;
1418 #endif
1419  Vector<T> Xi;
1420  X.refCol(i,Xi);
1421  Vector<bool> maski;
1422  mask.refCol(i,maski);
1423  Vector<int> ind;
1424  rM.refCol(i,ind);
1425  Vector<T> coeffs;
1426  vM.refCol(i,coeffs);
1427  Vector<T>& DtR=DtRT[numT];
1428 
1429  if (maski.allfalse()) continue;
1430  if (maski.alltrue()) {
1431  T normX = Xi.nrm2sq();
1432  D.multTrans(Xi,DtR);
1433  coreLARS2(DtR,G,GsT[numT],GaT[numT],invGsT[numT],uT[numT],coeffs,
1434  ind,workT[numT],normX,mode,constraint,pos);
1435  } else {
1436  D.copyMask(DmaskT[numT],maski);
1437  Xi.copyMask(XmaskT[numT],maski);
1438  T constraint_mask = mode == PENALTY || mode == L2ERROR ? constraint*XmaskT[numT].n()/Xi.n() : constraint;
1439  T normX = XmaskT[numT].nrm2sq();
1440  DmaskT[numT].multTrans(XmaskT[numT],DtR);
1441  GT[numT].setMatrices(DmaskT[numT],false);
1442  GT[numT].addDiag(MAX(lambda2,T(1e-10)));
1443  coreLARS2(DtR,GT[numT],
1444  GsT[numT],GaT[numT],invGsT[numT],uT[numT],coeffs,
1445  ind,workT[numT],normX,mode,constraint_mask,pos);
1446  DmaskT[numT].setm(D.m());
1447  DmaskT[numT].setn(D.n());
1448  XmaskT[numT].setn(X.m());
1449  }
1450  }
1451 
1452  delete[](GT);
1453  delete[](XmaskT);
1454  delete[](DmaskT);
1455  delete[](DtRT);
1456  delete[](uT);
1457  delete[](GsT);
1458  delete[](GaT);
1459  delete[](invGsT);
1460  delete[](workT);
1461 
1463  spalpha.convert(vM,rM,K);
1464 
1465 };
1466 
1467 template <typename T>
1468 void lasso2(const Matrix<T>& X, const Matrix<T>& D, SpMatrix<T>& spalpha,
1469  int L, const T constraint, const T lambda2, constraint_type mode, const bool pos,
1470  const int numThreads, Matrix<T>* path, int length_path) {
1471  ProdMatrix<T> G(D,X.n() > 10 && D.n() < 50000);
1472  ProdMatrix<T> DtX(D,X,false);
1473  G.addDiag(MAX(lambda2,1e-10));
1474  lasso2(X,G,DtX,spalpha,L,constraint,mode,pos,numThreads,path, length_path);
1475 }
1476 
1477 
1478 template <typename T>
1479 void lasso2(const Data<T>& X, const AbstractMatrix<T>& G, const AbstractMatrix<T>& DtX,
1480  SpMatrix<T>& spalpha,
1481  int L, const T constraint, constraint_type mode, const bool pos,
1482  const int numThreads, Matrix<T>* path, int length_path) {
1483  spalpha.clear();
1484  const int M = X.n();
1485  const int K = G.n();
1486  Matrix<T> vM;
1487  Matrix<int> rM;
1488  vM.resize(L,M);
1489  rM.resize(L,M);
1490 
1491  if (L <= 0) return;
1492  if (path) path->setZeros();
1493 
1494  int NUM_THREADS=init_omp(numThreads);
1495 
1496  Vector<T>* DtRT=new Vector<T>[NUM_THREADS];
1497  Vector<T>* uT=new Vector<T>[NUM_THREADS];
1498  Matrix<T>* GsT=new Matrix<T>[NUM_THREADS];
1499  Matrix<T>* GaT=new Matrix<T>[NUM_THREADS];
1500  Matrix<T>* invGsT=new Matrix<T>[NUM_THREADS];
1501  Matrix<T>* workT=new Matrix<T>[NUM_THREADS];
1502  for (int i = 0; i<NUM_THREADS; ++i) {
1503  DtRT[i].resize(K);
1504  uT[i].resize(K);
1505  uT[i].setZeros();
1506  GsT[i].resize(L,L);
1507  invGsT[i].resize(L,L);
1508  GaT[i].resize(K,L);
1509  workT[i].resize(K,3);
1510  workT[i].setZeros();
1511  }
1512  int i;
1513  Vector<T> norms;
1514  X.norm_2sq_cols(norms);
1515 #pragma omp parallel for private(i)
1516  for (i = 0; i< M; ++i) {
1517 #ifdef _OPENMP
1518  int numT=omp_get_thread_num();
1519 #else
1520  int numT=0;
1521 #endif
1522  // Vector<T> Xi;
1523  // X.refCol(i,Xi);
1524  // T normX = Xi.nrm2sq();
1525  T normX = norms[i];
1526 
1527  Vector<int> ind;
1528  rM.refCol(i,ind);
1529  Vector<T> coeffs;
1530  vM.refCol(i,coeffs);
1531 
1532  Vector<T>& DtR=DtRT[numT];
1533  DtX.copyCol(i,DtR);
1534  //D.multTrans(Xi,DtR);
1535  coreLARS2(DtR,G,GsT[numT],GaT[numT],invGsT[numT],
1536  uT[numT],coeffs,
1537  ind,workT[numT],normX,mode,constraint,pos,
1538  path && i==0 ? path->rawX() : NULL,length_path);
1539  }
1540 
1541  delete[](DtRT);
1542  delete[](uT);
1543  delete[](GsT);
1544  delete[](GaT);
1545  delete[](invGsT);
1546  delete[](workT);
1547 
1549  spalpha.convert(vM,rM,K);
1550 };
1551 
1552 
1553 
1555 template <typename T>
1556 void coreLARS2(Vector<T>& DtR, const AbstractMatrix<T>& G,
1557  Matrix<T>& Gs,
1558  Matrix<T>& Ga,
1559  Matrix<T>& invGs,
1560  Vector<T>& u,
1561  Vector<T>& coeffs,
1562  Vector<int>& ind,
1563  Matrix<T>& work,
1564  T& normX,
1565  const constraint_type mode,
1566  const T constraint,
1567  const bool pos,
1568  T* path, int length_path) {
1569  const int LL = Gs.n();
1570  const int K = G.n();
1571  const int L = MIN(LL,K);
1572  if (length_path <= 1) length_path=4*L;
1573 
1574  coeffs.setZeros();
1575  ind.set(-1);
1576 
1577  T* const pr_Gs = Gs.rawX();
1578  T* const pr_invGs = invGs.rawX();
1579  T* const pr_Ga = Ga.rawX();
1580  T* const pr_work = work.rawX();
1581  T* const pr_u = u.rawX();
1582  T* const pr_DtR = DtR.rawX();
1583  T* const pr_coeffs = coeffs.rawX();
1584  int* const pr_ind = ind.rawX();
1585 
1586  // Find the most correlated element
1587  int currentInd = pos ? DtR.max() : DtR.fmax();
1588  if (mode == PENALTY && abs(DtR[currentInd]) < constraint) return;
1589  if (mode == L2ERROR && normX < constraint) return;
1590  bool newAtom=true;
1591 
1592  int i;
1593  int iter=0;
1594  T thrs = 0;
1595  for (i = 0; i<L; ++i) {
1596  ++iter;
1597  if (newAtom) {
1598  pr_ind[i]=currentInd;
1599  // cerr << "Add " << currentInd << endl;
1600  G.extract_rawCol(pr_ind[i],pr_Ga+i*K);
1601  for (int j = 0; j<=i; ++j)
1602  pr_Gs[i*LL+j]=pr_Ga[i*K+pr_ind[j]];
1603 
1604  // Update inverse of Gs
1605  if (i == 0) {
1606  pr_invGs[0]=T(1.0)/pr_Gs[0];
1607  } else {
1608  cblas_symv<T>(CblasColMajor,CblasUpper,i,T(1.0),
1609  pr_invGs,LL,pr_Gs+i*LL,1,T(0.0),pr_u,1);
1610  const T schur =
1611  T(1.0)/(pr_Gs[i*LL+i]-cblas_dot<T>(i,pr_u,1,pr_Gs+i*LL,1));
1612  pr_invGs[i*LL+i]=schur;
1613  cblas_copy<T>(i,pr_u,1,pr_invGs+i*LL,1);
1614  cblas_scal<T>(i,-schur,pr_invGs+i*LL,1);
1615  cblas_syr<T>(CblasColMajor,CblasUpper,i,schur,pr_u,1,
1616  pr_invGs,LL);
1617  }
1618  }
1619 
1620  // Compute the path direction
1621  for (int j = 0; j<=i; ++j)
1622  pr_work[j]= pr_DtR[pr_ind[j]] > 0 ? T(1.0) : T(-1.0);
1623  cblas_symv<T>(CblasColMajor,CblasUpper,i+1,T(1.0),pr_invGs,LL,
1624  pr_work,1,T(0.0),pr_u,1);
1625 
1626  // Compute the step on the path
1627  T step_max = INFINITY;
1628  int first_zero = -1;
1629  for (int j = 0; j<=i; ++j) {
1630  T ratio = -pr_coeffs[j]/pr_u[j];
1631  if (ratio > 0 && ratio <= step_max) {
1632  step_max=ratio;
1633  first_zero=j;
1634  }
1635  }
1636  // PRINT_F(step_max)
1637 
1638  T current_correlation = abs<T>(pr_DtR[pr_ind[0]]);
1639  cblas_gemv<T>(CblasColMajor,CblasNoTrans,K,i+1,T(1.0),pr_Ga,
1640  K,pr_u,1,T(0.0),pr_work+2*K,1);
1641  cblas_copy<T>(K,pr_work+2*K,1,pr_work+K,1);
1642  cblas_copy<T>(K,pr_work+2*K,1,pr_work,1);
1643 
1644  for (int j = 0; j<=i; ++j) {
1645  pr_work[pr_ind[j]]=INFINITY;
1646  pr_work[pr_ind[j]+K]=INFINITY;
1647  }
1648  for (int j = 0; j<K; ++j) {
1649  pr_work[j] = ((pr_work[j] < INFINITY) && (pr_work[j] > T(-1.0))) ? (pr_DtR[j]+current_correlation)/(T(1.0)+pr_work[j]) : INFINITY;
1650  }
1651  // work.print("work");
1652  for (int j = 0; j<K; ++j) {
1653  pr_work[j+K] = ((pr_work[j+K] < INFINITY) && (pr_work[j+K] < T(1.0))) ? (current_correlation-pr_DtR[j])/(T(1.0)-pr_work[j+K]) : INFINITY;
1654  }
1655  // work.print("work");
1656 
1657  if (pos) {
1658  for (int j = 0; j<K; ++j) {
1659  pr_work[j]=INFINITY;
1660  }
1661  }
1662  // work.print("work");
1663  // coeffs.print("coeffs");
1664  int index = cblas_iamin<T>(2*K,pr_work,1);
1665  T step = pr_work[index];
1666 
1667  // Choose next element
1668  currentInd = index % K;
1669 
1670  // compute the coefficients of the polynome representing normX^2
1671  T coeff1 = 0;
1672  for (int j = 0; j<=i; ++j)
1673  coeff1 += pr_DtR[pr_ind[j]] > 0 ? pr_u[j] : -pr_u[j];
1674  T coeff2 = 0;
1675  for (int j = 0; j<=i; ++j)
1676  coeff2 += pr_DtR[pr_ind[j]]*pr_u[j];
1677  T coeff3 = normX-constraint;
1678 
1679 
1680  T step_max2;
1681  if (mode == PENALTY) {
1682  step_max2 = current_correlation-constraint;
1683  } else if (mode == L2ERROR) {
1685  const T delta = coeff2*coeff2-coeff1*coeff3;
1686  step_max2 = delta < 0 ? INFINITY : (coeff2-sqrt(delta))/coeff1;
1687  step_max2 = MIN(current_correlation,step_max2);
1688  } else {
1690  step_max2 = coeff1 < 0 ? INFINITY : (constraint-thrs)/coeff1;
1691  step_max2 = MIN(current_correlation,step_max2);
1692  }
1693  step = MIN(MIN(step,step_max2),step_max);
1694  if (step == INFINITY) break; // stop the path
1695 
1696  // Update coefficients
1697  cblas_axpy<T>(i+1,step,pr_u,1,pr_coeffs,1);
1698 
1699  if (pos) {
1700  for (int j = 0; j<i+1; ++j)
1701  if (pr_coeffs[j] < 0) pr_coeffs[j]=0;
1702  }
1703 
1704  // Update correlations
1705  cblas_axpy<T>(K,-step,pr_work+2*K,1,pr_DtR,1);
1706 
1707  // Update normX
1708  normX += coeff1*step*step-2*coeff2*step;
1709 
1710  // Update norm1
1711  thrs += step*coeff1;
1712 
1713  if (path) {
1714  for (int k = 0; k<=i; ++k)
1715  path[iter*K+ind[k]]=pr_coeffs[k];
1716  }
1717 
1718  // Choose next action
1719 
1720  if (step == step_max) {
1721  // cerr << "Remove " << pr_ind[first_zero] << endl;
1724  for (int j = first_zero; j<i; ++j) {
1725  cblas_copy<T>(K,pr_Ga+(j+1)*K,1,pr_Ga+j*K,1);
1726  pr_ind[j]=pr_ind[j+1];
1727  pr_coeffs[j]=pr_coeffs[j+1];
1728  }
1729  pr_ind[i]=-1;
1730  pr_coeffs[i]=0;
1731  for (int j = first_zero; j<i; ++j) {
1732  cblas_copy<T>(first_zero,pr_Gs+(j+1)*LL,1,pr_Gs+j*LL,1);
1733  cblas_copy<T>(i-first_zero,pr_Gs+(j+1)*LL+first_zero+1,1,
1734  pr_Gs+j*LL+first_zero,1);
1735  }
1736  const T schur = pr_invGs[first_zero*LL+first_zero];
1737  cblas_copy<T>(first_zero,pr_invGs+first_zero*LL,1,pr_u,1);
1738  cblas_copy<T>(i-first_zero,pr_invGs+(first_zero+1)*LL+first_zero,LL,
1739  pr_u+first_zero,1);
1740  for (int j = first_zero; j<i; ++j) {
1741  cblas_copy<T>(first_zero,pr_invGs+(j+1)*LL,1,pr_invGs+j*LL,1);
1742  cblas_copy<T>(i-first_zero,pr_invGs+(j+1)*LL+first_zero+1,1,
1743  pr_invGs+j*LL+first_zero,1);
1744  }
1745  cblas_syr<T>(CblasColMajor,CblasUpper,i,T(-1.0)/schur,
1746  pr_u,1,pr_invGs,LL);
1747  newAtom=false;
1748  i=i-2;
1749  } else {
1750  newAtom=true;
1751  }
1752  if ((iter >= length_path-1) || abs(step) < 1e-15 ||
1753  step == step_max2 || (normX < 1e-15) ||
1754  (i == (L-1)) ||
1755  (mode == L2ERROR && normX - constraint < 1e-15) ||
1756  (mode == L1COEFFS && (constraint-thrs < 1e-15))) {
1757  break;
1758  }
1759  }
1760 }
1761 
1763 template <typename T>
1764 void coreLARS2W(Vector<T>& DtR, AbstractMatrix<T>& G,
1765  Matrix<T>& Gs,
1766  Matrix<T>& Ga,
1767  Matrix<T>& invGs,
1768  Vector<T>& u,
1769  Vector<T>& coeffs,
1770  const Vector<T>& weights,
1771  Vector<int>& ind,
1772  Matrix<T>& work,
1773  T& normX,
1774  const constraint_type mode,
1775  const T constraint,
1776  const bool pos) {
1777  const int LL = Gs.n();
1778  const int K = G.n();
1779  const int L = MIN(LL,K);
1780  coeffs.setZeros();
1781  ind.set(-1);
1782 
1783  T* const pr_Gs = Gs.rawX();
1784  T* const pr_invGs = invGs.rawX();
1785  T* const pr_Ga = Ga.rawX();
1786  // T* const pr_G = G.rawX();
1787  T* const pr_work = work.rawX();
1788  T* const pr_u = u.rawX();
1789  T* const pr_DtR = DtR.rawX();
1790  T* const pr_coeffs = coeffs.rawX();
1791  T* const pr_weights = weights.rawX();
1792  int* const pr_ind = ind.rawX();
1793 
1794  DtR.div(weights);
1795 
1796  // Find the most correlated element
1797  int currentInd = pos ? DtR.max() : DtR.fmax();
1798  if (mode == PENALTY && abs(DtR[currentInd]) < constraint) return;
1799  if (mode == L2ERROR && normX < constraint) return;
1800  bool newAtom=true;
1801 
1802  int i;
1803  int iter=0;
1804  T thrs = 0;
1805  for (i = 0; i<L; ++i) {
1806  ++iter;
1807  if (newAtom) {
1808  pr_ind[i]=currentInd;
1809  // Update upper part of Gs and Ga
1810  G.extract_rawCol(pr_ind[i],pr_Ga+i*K);
1811  for (int j = 0; j<=i; ++j)
1812  pr_Gs[i*LL+j]=pr_Ga[i*K+pr_ind[j]];
1813 
1814  // Update inverse of Gs
1815  if (i == 0) {
1816  pr_invGs[0]=T(1.0)/pr_Gs[0];
1817  } else {
1818  cblas_symv<T>(CblasColMajor,CblasUpper,i,T(1.0),
1819  pr_invGs,LL,pr_Gs+i*LL,1,T(0.0),pr_u,1);
1820  const T schur =
1821  T(1.0)/(pr_Gs[i*LL+i]-cblas_dot<T>(i,pr_u,1,pr_Gs+i*LL,1));
1822  pr_invGs[i*LL+i]=schur;
1823  cblas_copy<T>(i,pr_u,1,pr_invGs+i*LL,1);
1824  cblas_scal<T>(i,-schur,pr_invGs+i*LL,1);
1825  cblas_syr<T>(CblasColMajor,CblasUpper,i,schur,pr_u,1,
1826  pr_invGs,LL);
1827  }
1828  }
1829 
1830  // Compute the path direction
1831  for (int j = 0; j<=i; ++j)
1832  pr_work[j]= pr_DtR[pr_ind[j]] > 0 ? weights[pr_ind[j]] : -weights[pr_ind[j]];
1833  cblas_symv<T>(CblasColMajor,CblasUpper,i+1,T(1.0),pr_invGs,LL,
1834  pr_work,1,T(0.0),pr_u,1);
1835 
1836  // Compute the step on the path
1837  T step_max = INFINITY;
1838  int first_zero = -1;
1839  for (int j = 0; j<=i; ++j) {
1840  T ratio = -pr_coeffs[j]/pr_u[j];
1841  if (ratio > 0 && ratio <= step_max) {
1842  step_max=ratio;
1843  first_zero=j;
1844  }
1845  }
1846 
1847  T current_correlation = abs<T>(pr_DtR[pr_ind[0]]);
1848  cblas_gemv<T>(CblasColMajor,CblasNoTrans,K,i+1,T(1.0),pr_Ga,
1849  K,pr_u,1,T(0.0),pr_work+2*K,1);
1850  vDiv<T>(K,pr_work+2*K,pr_weights,pr_work+2*K);
1851  cblas_copy<T>(K,pr_work+2*K,1,pr_work+K,1);
1852  cblas_copy<T>(K,pr_work+2*K,1,pr_work,1);
1853 
1854  for (int j = 0; j<=i; ++j) {
1855  pr_work[pr_ind[j]]=INFINITY;
1856  pr_work[pr_ind[j]+K]=INFINITY;
1857  }
1858  for (int j = 0; j<K; ++j) {
1859  pr_work[j] = ((pr_work[j] < INFINITY) && (pr_work[j] > T(-1.0))) ? (pr_DtR[j]+current_correlation)/(T(1.0)+pr_work[j]) : INFINITY;
1860  }
1861  for (int j = 0; j<K; ++j) {
1862  pr_work[j+K] = ((pr_work[j+K] < INFINITY) && (pr_work[j+K] < T(1.0))) ? (current_correlation-pr_DtR[j])/(T(1.0)-pr_work[j+K]) : INFINITY;
1863  }
1864 
1865  if (pos) {
1866  for (int j = 0; j<K; ++j) {
1867  pr_work[j]=INFINITY;
1868  }
1869  }
1870  int index = cblas_iamin<T>(2*K,pr_work,1);
1871  T step = pr_work[index];
1872  // Choose next element
1873  currentInd = index % K;
1874 
1875  // compute the coefficients of the polynome representing normX^2
1876  T coeff1 = 0;
1877  for (int j = 0; j<=i; ++j)
1878  coeff1 += pr_DtR[pr_ind[j]] > 0 ? pr_weights[pr_ind[j]]*pr_u[j] :
1879  -pr_weights[pr_ind[j]]*pr_u[j];
1880  T coeff2 = 0;
1881  for (int j = 0; j<=i; ++j)
1882  coeff2 += pr_DtR[pr_ind[j]]*pr_u[j]*pr_weights[pr_ind[j]];
1883  T coeff3 = normX-constraint;
1884 
1885  T step_max2;
1886  if (mode == PENALTY) {
1887  step_max2 = current_correlation-constraint;
1888  } else if (mode == L2ERROR) {
1890  const T delta = coeff2*coeff2-coeff1*coeff3;
1891  step_max2 = delta < 0 ? INFINITY : (coeff2-sqrt(delta))/coeff1;
1892  } else {
1894  step_max2 = coeff1 < 0 ? INFINITY : (constraint-thrs)/coeff1;
1895  }
1896  step = MIN(MIN(step,step_max2),step_max);
1897 
1898  if (step == INFINITY) break; // stop the path
1899 
1900  // Update coefficients
1901  cblas_axpy<T>(i+1,step,pr_u,1,pr_coeffs,1);
1902 
1903  // Update correlations
1904  cblas_axpy<T>(K,-step,pr_work+2*K,1,pr_DtR,1);
1905 
1906  // Update normX
1907  normX += coeff1*step*step-2*coeff2*step;
1908 
1909  // Update norm1
1910  thrs += step*coeff1;
1911 
1912  if (step == step_max) {
1915  for (int j = first_zero; j<i; ++j) {
1916  cblas_copy<T>(K,pr_Ga+(j+1)*K,1,pr_Ga+j*K,1);
1917  pr_ind[j]=pr_ind[j+1];
1918  pr_coeffs[j]=pr_coeffs[j+1];
1919  }
1920  pr_ind[i]=-1;
1921  pr_coeffs[i]=0;
1922  for (int j = first_zero; j<i; ++j) {
1923  cblas_copy<T>(first_zero,pr_Gs+(j+1)*LL,1,pr_Gs+j*LL,1);
1924  cblas_copy<T>(i-first_zero,pr_Gs+(j+1)*LL+first_zero+1,1,
1925  pr_Gs+j*LL+first_zero,1);
1926  }
1927  const T schur = pr_invGs[first_zero*LL+first_zero];
1928  cblas_copy<T>(first_zero,pr_invGs+first_zero*LL,1,pr_u,1);
1929  cblas_copy<T>(i-first_zero,pr_invGs+(first_zero+1)*LL+first_zero,LL,
1930  pr_u+first_zero,1);
1931  for (int j = first_zero; j<i; ++j) {
1932  cblas_copy<T>(first_zero,pr_invGs+(j+1)*LL,1,pr_invGs+j*LL,1);
1933  cblas_copy<T>(i-first_zero,pr_invGs+(j+1)*LL+first_zero+1,1,
1934  pr_invGs+j*LL+first_zero,1);
1935  }
1936  cblas_syr<T>(CblasColMajor,CblasUpper,i,T(-1.0)/schur,
1937  pr_u,1,pr_invGs,LL);
1938  newAtom=false;
1939  i=i-2;
1940  } else {
1941  newAtom=true;
1942  }
1943  // Choose next action
1944  if (iter > 4*L || abs(step) < 1e-10 ||
1945  step == step_max2 || (normX < 1e-10) ||
1946  (i == (L-1)) ||
1947  (mode == L2ERROR && normX - constraint < 1e-10) ||
1948  (mode == L1COEFFS && (constraint-thrs < 1e-10))) {
1949  break;
1950  }
1951  }
1952 }
1953 
1954 
1955 
1956 /* ************************
1957  * Iterative thresholding
1958  * ************************/
1959 
1965 template <typename T>
1966 void ist(const Matrix<T>& X, const Matrix<T>& D,
1967  SpMatrix<T>& spalpha, T lambda, constraint_type mode,
1968  const int itermax,
1969  const T tol,
1970  const int numThreads) {
1971  Matrix<T> alpha;
1972  spalpha.toFull(alpha);
1973  spalpha.clear();
1974  ist(X,D,alpha,lambda,mode,itermax,tol,numThreads);
1975  alpha.toSparse(spalpha);
1976 }
1977 
1978 template <typename T>
1979 void ist(const Matrix<T>& X, const Matrix<T>& D,
1980  Matrix<T>& alpha, T lambda, constraint_type mode,
1981  const int itermax,
1982  const T tol, const int numThreads) {
1983 
1984  if (mode == L1COEFFS) {
1985  std::cerr << "Mode not implemented" << std::endl;
1986  return;
1987  }
1988 
1989  int K=D.n();
1990  int M=X.n();
1991  alpha.resize(K,M);
1992  if (!D.isNormalized()) {
1993  cerr << "Current implementation of IST does not support non-normalized dictionaries" << endl;
1994  return;
1995  }
1996 
1998  //CachedProdMatrix<T> G(D, K < 20000 && M*K/10 > K);
1999  //ProdMatrix<T> G(D, K < 20000 && M*K/10 > K);
2000  Matrix<T> G;
2001  D.XtX(G);
2002  // for (int i = 0; i<K; ++i) G[i*K+i] += 1e-6;
2003  G.addDiag(1e-12);
2004  ProdMatrix<T> DtX(D,X,false);
2005 
2006  int NUM_THREADS=init_omp(numThreads);
2007 
2008  Vector<T>* DtRT= new Vector<T>[NUM_THREADS];
2009  SpVector<T>* spAlphaT= new SpVector<T>[NUM_THREADS];
2010  for (int i = 0; i<NUM_THREADS; ++i) {
2011  DtRT[i].resize(K);
2012  spAlphaT[i].resize(K);
2013  };
2014 
2015  int i;
2016 #pragma omp parallel for private(i)
2017  for (i = 0; i< M; ++i) {
2018 #ifdef _OPENMP
2019  int numT=omp_get_thread_num();
2020 #else
2021  int numT=0;
2022 #endif
2023  Vector<T> coeffs;
2024  alpha.refCol(i,coeffs);
2025  Vector<T>& DtR=DtRT[numT];
2026  SpVector<T>& spAlpha=spAlphaT[numT];
2027  T norm1 = coeffs.asum();
2028  // Compute DtR
2029  DtX.copyCol(i,DtR);
2030  Vector<T> Xi;
2031  X.refCol(i,Xi);
2032  T normX2 = Xi.nrm2sq();
2033 
2034  if (norm1 > EPSILON) {
2035  coeffs.toSparse(spAlpha);
2036  G.mult(spAlpha,DtR,-1.0,1.0);
2037  }
2038 
2039  if (mode == PENALTY) {
2040  coreIST(G,DtR,coeffs,lambda,itermax,tol);
2041  } else {
2042  coreISTconstrained(G,DtR,coeffs,normX2,lambda,itermax,tol);
2043  }
2044  }
2045 
2046  delete[](DtRT);
2047  delete[](spAlphaT);
2048 
2049 }
2050 
2051 template <typename T>
2052 inline void coreIST(const AbstractMatrix<T>& G, Vector<T>& DtRv, Vector<T>& coeffsv,
2053  const T thrs, const int itermax,
2054  const T tol) {
2055 
2056  const int K = G.n();
2057  T* const coeffs = coeffsv.rawX();
2058  T* const DtR = DtRv.rawX();
2059  // T* const prG = G.rawX();
2060 
2061  const T lambda_init=thrs;
2062  T maxDtR = DtRv.fmaxval();
2063  T norm1=coeffsv.asum();
2064  T lambda=lambda_init;
2065  vAdd(K,DtR,coeffs,DtR);
2066 
2067  for (int iter=0; iter < itermax; ++iter) {
2068  for (int j = 0; j <K; ++j) {
2069  if (DtR[j] > lambda) {
2070  T diff=coeffs[j];
2071  coeffs[j]=DtR[j]-lambda;
2072  diff-=coeffs[j];
2073  DtR[j]-=diff;
2074  G.add_rawCol(j,DtR,diff);
2075  //cblas_axpy(K,diff,prG+j*K,1,DtR,1);
2076  } else if (DtR[j] < -lambda) {
2077  T diff=coeffs[j];
2078  coeffs[j]=DtR[j]+lambda;
2079  diff-=coeffs[j];
2080  DtR[j]-=diff;
2081  G.add_rawCol(j,DtR,diff);
2082  //cblas_axpy(K,diff,prG+j*K,1,DtR,1);
2083  } else if (coeffs[j]) {
2084  T diff=coeffs[j];
2085  coeffs[j]=T();
2086  DtR[j]-=diff;
2087  G.add_rawCol(j,DtR,diff);
2088  //cblas_axpy(K,diff,prG+j*K,1,DtR,1);
2089  }
2090  }
2091  if (iter % 5 == 1) {
2092  vSub(K,DtR,coeffs,DtR);
2093  maxDtR = DtRv.fmaxval();
2094  norm1 =T();
2095  T DtRa = T();
2096  for (int j = 0; j<K; ++j) {
2097  if (coeffs[j]) {
2098  norm1 += abs(coeffs[j]);
2099  DtRa += DtR[j]*coeffs[j];
2100  }
2101  }
2102  vAdd(K,DtR,coeffs,DtR);
2103  const T kappa = -DtRa+norm1*maxDtR;
2104  if (abs(lambda - maxDtR) < tol && kappa <= tol)
2105  break;
2106  }
2107  }
2108 }
2109 
2110 
2112 template <typename T>
2113 void coreISTconstrained(const AbstractMatrix<T>& G, Vector<T>& DtRv, Vector<T>&
2114  coeffsv, const T normX2, const T eps, const int itermax, const T tol) {
2115  const int K = G.n();
2116  T* const coeffs = coeffsv.rawX();
2117  T* const DtR = DtRv.rawX();
2118  // T* const prG = G.rawX();
2119  T err = normX2;
2120 
2121  T norm1 = coeffsv.asum();
2122  if (!norm1 && err <= eps) return;
2123  T current_tol = 10.0*tol;
2124  T maxDtR = DtRv.fmaxval();
2125  T lambda = maxDtR;
2126  T lambdasq= lambda*lambda;
2127  if (!norm1) {
2128  lambdasq *= eps/err;
2129  lambda=sqrt(lambdasq);
2130  }
2131 
2132  Vector<int> indices(K);
2133  indices.set(-1);
2134  int* const pr_indices=indices.rawX();
2135  int count;
2136 
2137  for (int iter=0; iter < itermax; ++iter) {
2138 
2139  count=0;
2140  T old_err = err;
2141  for (int j = 0; j <K; ++j) {
2142 
2143  // Soft-thresholding
2144  T old_coeff = coeffs[j];
2145  T diff = DtR[j]+old_coeff;
2146  if (diff > lambda) {
2147  coeffs[j] = diff - lambda;
2148  err+=lambdasq-DtR[j]*DtR[j];
2149  pr_indices[count++]=j;
2150  } else if (diff < - lambda) {
2151  coeffs[j] = diff + lambda;
2152  err+=lambdasq-DtR[j]*DtR[j];
2153  pr_indices[count++]=j;
2154  } else {
2155  coeffs[j]=T();
2156  if (old_coeff) {
2157  err+=diff*diff-DtR[j]*DtR[j];
2158  }
2159  }
2160  // Update DtR
2161  diff = old_coeff-coeffs[j];
2162  if (diff) {
2163  G.add_rawCol(j,DtR,diff);
2164  //cblas_axpy<T>(K,old_coeff-coeffs[j],prG+j*K,1,DtR,1);
2165  }
2166  }
2167 
2168  maxDtR = DtRv.fmaxval();
2169  norm1 =T();
2170  T DtRa = T();
2171  for (int j = 0; j<count; ++j) {
2172  const int ind = pr_indices[j];
2173  norm1 += abs(coeffs[ind]);
2174  DtRa += DtR[ind]*coeffs[ind];
2175  }
2176  if (norm1-DtRa/maxDtR <= current_tol) {
2177  const bool change = ((old_err > eps) && err < eps+current_tol) ||
2178  (old_err < eps && err > eps-current_tol);
2179  if (change) {
2180  if (current_tol == tol) {
2181  break;
2182  } else {
2183  current_tol = MAX(current_tol*0.5,tol);
2184  }
2185  }
2186  lambdasq *= eps/err;
2187  lambda=sqrt(lambdasq);
2188  }
2189  }
2190 };
2191 
2192 
2193 
2195 template <typename T>
2196 void ist_groupLasso(const Matrix<T>* XT, const Matrix<T>& D,
2197  Matrix<T>* alphaT, const int Ngroups,
2198  const T lambda, const constraint_type mode,
2199  const int itermax,
2200  const T tol, const int numThreads) {
2201  int K=D.n();
2202  int n = D.m();
2203 
2204  if (!D.isNormalized()) {
2205  cerr << "Current implementation of block coordinate descent does not support non-normalized dictionaries" << endl;
2206  return;
2207  }
2208 
2209  if (mode == L1COEFFS) {
2210  std::cerr << "Mode not implemented" << std::endl;
2211  return;
2212  }
2213 
2214 
2216  Matrix<T> G;
2217  D.XtX(G);
2218 
2219  int NUM_THREADS=init_omp(numThreads);
2220 
2221  Matrix<T>* RtDT = new Matrix<T>[NUM_THREADS];
2222  Matrix<T>* alphatT = new Matrix<T>[NUM_THREADS];
2223 
2224  int i;
2225 #pragma omp parallel for private(i)
2226  for (i = 0; i< Ngroups; ++i) {
2227 #ifdef _OPENMP
2228  int numT=omp_get_thread_num();
2229 #else
2230  int numT=0;
2231 #endif
2232  const Matrix<T>& X = XT[i];
2233  int M = X.n();
2234  Matrix<T>& alphat = alphatT[numT];
2235  alphaT[i].transpose(alphat);
2236  Matrix<T>& RtD = RtDT[numT];
2237  X.mult(D,RtD,true,false);
2238 
2239 
2240  Vector<T> col, col2;
2241  T norm1 = alphat.asum();
2242  T normX2;
2243 
2244  if (!norm1) {
2245  Vector<T> DtR_mean(K);
2246  Vector<T> coeffs_mean(K);
2247  coeffs_mean.setZeros();
2248  RtD.meanRow(DtR_mean);
2249  coeffs_mean.setZeros();
2250  if (mode == PENALTY) {
2251  coreIST(G,DtR_mean,coeffs_mean,lambda/T(2.0),itermax,tol);
2252  } else {
2253  Vector<T> meanVec(n);
2254  X.meanCol(meanVec);
2255  normX2=meanVec.nrm2sq();
2256  coreISTconstrained(G,DtR_mean,coeffs_mean,normX2,
2257  lambda,itermax,tol);
2258  SpVector<T> spalpha(K);
2259  normX2-=computeError(normX2,G,DtR_mean,coeffs_mean,spalpha);
2260  normX2=X.normFsq()-M*normX2;
2261  }
2262  alphat.fillRow(coeffs_mean);
2263  }
2264 
2265  if (M > 1) {
2266  for (int j = 0; j<K; ++j) {
2267  alphat.refCol(j,col);
2268  const T nrm=col.nrm2sq();
2269  if (nrm) {
2270  G.refCol(j,col2);
2271  RtD.rank1Update(col,col2,T(-1.0));
2272  }
2273  }
2274 
2275  if (mode == PENALTY) {
2276  coreGroupIST(G,RtD,alphat,sqr<T>(M)*lambda/T(2.0),itermax,sqr<T>(M)*tol);
2277  } else {
2278  coreGroupISTConstrained(G,RtD,alphat,normX2,M*lambda,itermax,sqr<T>(M)*tol);
2279  }
2280  }
2281  alphat.transpose(alphaT[i]);
2282  }
2283 
2284  delete[](RtDT);
2285  delete[](alphatT);
2286 };
2287 
2288 
2289 template <typename T>
2290 void coreGroupIST(const Matrix<T>& G, Matrix<T>& RtDm,
2291  Matrix<T>& coeffsm,
2292  const T thrs,
2293  const int itermax,
2294  const T tol) {
2295  const int K = G.n();
2296  const int M = RtDm.m();
2297  T* const prG = G.rawX();
2298  T* const RtD = RtDm.rawX();
2299  T* const coeffs = coeffsm.rawX();
2300 
2301  const T lambda_init=thrs;
2302  T lambda=lambda_init;
2303 
2304  Vector<T> old_coeffv(M);
2305  T* const old_coeff = old_coeffv.rawX();
2306  Vector<T> normsv(K);
2307  T* const norms = normsv.rawX();
2308  coeffsm.norm_2_cols(normsv);
2309  Vector<T> normRtDv(K);
2310 
2311  Vector<int> activatev(K);
2312  activatev.set(3);
2313  int* const activate=activatev.rawX();
2314 
2315  for (int iter=0; iter < itermax; ++iter) {
2316  for (int j = 0; j <K; ++j) {
2317  if (activate[j] >= 0) {
2318  if (norms[j]) {
2319  cblas_copy(M,coeffs+j*M,1,old_coeff,1);
2320  vAdd(M,coeffs+j*M,RtD+j*M,coeffs+j*M);
2321  const T nrm = cblas_nrm2(M,coeffs+j*M,1);
2322  if (nrm > lambda) {
2323  norms[j]=nrm-lambda;
2324  cblas_scal(M,norms[j]/nrm,coeffs+j*M,1);
2325  vSub(M,old_coeff,coeffs+j*M,old_coeff);
2326  cblas_ger(CblasColMajor,M,K,T(1.0),old_coeff,1,prG+j*K,1,RtD,M);
2327  activate[j]=5;
2328  } else {
2329  memset(coeffs+j*M,0,M*sizeof(T));
2330  norms[j]=T();
2331  cblas_ger(CblasColMajor,M,K,T(1.0),old_coeff,1,prG+j*K,1,RtD,M);
2332  --activate[j];
2333  }
2334  } else {
2335  cblas_copy(M,RtD+j*M,1,old_coeff,1);
2336  const T nrm = cblas_nrm2(M,old_coeff,1);
2337  if (nrm > lambda) {
2338  norms[j]=nrm-lambda;
2339  cblas_copy(M,old_coeff,1,coeffs+j*M,1);
2340  cblas_scal(M,norms[j]/nrm,coeffs+j*M,1);
2341  cblas_ger(CblasColMajor,M,K,T(-1.0),coeffs+j*M,1,prG+j*K,1,RtD,M);
2342  activate[j]=5;
2343  } else {
2344  activate[j] = (activate[j] == 0) ? -10 : activate[j]-1;
2345  }
2346  }
2347  } else {
2348  ++activate[j];
2349  }
2350  }
2351 
2352  if (iter % 5 == 4) {
2353  T norm1=normsv.asum();
2354  RtDm.norm_2sq_cols(normRtDv);
2355  T maxDtR = sqr(normRtDv.maxval());
2356  T DtRa=T();
2357  for (int j = 0; j<K; ++j) {
2358  if (norms[j]) {
2359  DtRa += cblas_dot(M,coeffs+j*M,1,RtD+j*M,1);
2360  }
2361  }
2362  if ((maxDtR - lambda) < (tol*maxDtR/norm1) && norm1-DtRa/maxDtR < tol) break;
2363  }
2364  }
2365 };
2366 
2367 
2369 template <typename T>
2370 void coreGroupISTConstrained(const Matrix<T>& G, Matrix<T>& RtDm,
2371  Matrix<T>& coeffsm, const T normR,
2372  const T eps,
2373  const int itermax,
2374  const T tol) {
2375  const int K = G.n();
2376  const int M = RtDm.m();
2377  T* const prG = G.rawX();
2378  T* const RtD = RtDm.rawX();
2379  T* const coeffs = coeffsm.rawX();
2380 
2381  T err = normR;
2382 
2383  Vector<T> old_coeffv(M);
2384  T* const old_coeff = old_coeffv.rawX();
2385  Vector<T> normsv(K);
2386  T* const norms = normsv.rawX();
2387  coeffsm.norm_2_cols(normsv);
2388  Vector<T> normRtDv(K);
2389  RtDm.norm_2sq_cols(normRtDv);
2390 
2391  Vector<int> activatev(K);
2392  activatev.set(3);
2393  int* const activate=activatev.rawX();
2394 
2395  T norm1 = normsv.sum();
2396  if (!norm1 && err <= eps) return;
2397  T current_tol = 10.0*tol;
2398 
2399  T maxDtR = sqr(normRtDv.maxval());
2400  T lambda = maxDtR;
2401  T lambdasq= lambda*lambda;
2402 
2403  if (!norm1) {
2404  lambdasq *= eps/err;
2405  lambda=sqrt(lambdasq);
2406  }
2407 
2408  for (int iter=0; iter < itermax; ++iter) {
2409 
2410  T old_err = err;
2411  for (int j = 0; j <K; ++j) {
2412  if (activate[j] >= 0) {
2413  if (norms[j]) {
2414  cblas_copy(M,coeffs+j*M,1,old_coeff,1);
2415  vAdd(M,coeffs+j*M,RtD+j*M,coeffs+j*M);
2416  const T nrm = cblas_nrm2(M,coeffs+j*M,1);
2417  if (nrm > lambda) {
2418  norms[j]=nrm-lambda;
2419  cblas_scal(M,norms[j]/nrm,coeffs+j*M,1);
2420  vSub(M,old_coeff,coeffs+j*M,old_coeff);
2421  err += cblas_dot(M,old_coeff,1,old_coeff,1)
2422  +2*cblas_dot(M,old_coeff,1,RtD+j*M,1);
2423  cblas_ger(CblasColMajor,M,K,T(1.0),old_coeff,1,prG+j*K,1,RtD,M);
2424  activate[j]=3;
2425  } else {
2426  memset(coeffs+j*M,0,M*sizeof(T));
2427  norms[j]=T();
2428  err += cblas_dot(M,old_coeff,1,old_coeff,1)
2429  +2*cblas_dot(M,old_coeff,1,RtD+j*M,1);
2430  cblas_ger(CblasColMajor,M,K,T(1.0),old_coeff,1,prG+j*K,1,RtD,M);
2431  --activate[j];
2432  }
2433  } else {
2434  cblas_copy(M,RtD+j*M,1,old_coeff,1);
2435  const T nrm = cblas_nrm2(M,old_coeff,1);
2436  if (nrm > lambda) {
2437  norms[j]=nrm-lambda;
2438  cblas_copy(M,old_coeff,1,coeffs+j*M,1);
2439  cblas_scal(M,norms[j]/nrm,coeffs+j*M,1);
2440  err += cblas_dot(M,coeffs+j*M,1,coeffs+j*M,1)
2441  -2*cblas_dot(M,coeffs+j*M,1,RtD+j*M,1);
2442  cblas_ger(CblasColMajor,M,K,T(-1.0),coeffs+j*M,1,prG+j*K,1,RtD,M);
2443  activate[j]=3;
2444  } else {
2445  activate[j] = (activate[j] == 0) ? -3 : activate[j]-1;
2446  }
2447  }
2448  } else {
2449  ++activate[j];
2450  }
2451  }
2452 
2453  norm1 = normsv.sum();
2454  RtDm.norm_2sq_cols(normRtDv);
2455  maxDtR = sqr(normRtDv.maxval());
2456  T DtRa=T();
2457  for (int j = 0; j<K; ++j) {
2458  if (norms[j]) {
2459  DtRa += cblas_dot(M,coeffs+j*M,1,RtD+j*M,1);
2460  }
2461  }
2462  if (norm1-DtRa/maxDtR <= current_tol) {
2463  const T tol_bis=current_tol*maxDtR;
2464  const bool change = ((old_err > eps) && err < eps+tol_bis) ||
2465  (old_err < eps && err > eps-tol_bis);
2466  if (change) {
2467  if (current_tol == tol) {
2468  break;
2469  } else {
2470  current_tol = MAX(current_tol*0.5,tol);
2471  }
2472  }
2473  lambdasq *= eps/err;
2474  lambda=sqrt(lambdasq);
2475  }
2476  }
2477 };
2478 
2480 template <typename T>
2481 T computeError(const T normX2,const Vector<T>& norms,
2482  const Matrix<T>& G,const Matrix<T>& RtD,const Matrix<T>& alphat) {
2483  T err2 = normX2;
2484  Vector<T> col,col2;
2485  for (int j = 0; j<G.n(); ++j) {
2486  if (norms[j] > EPSILON) {
2487  alphat.refCol(j,col);
2488  RtD.refCol(j,col2);
2489  err2 -= 2*col.dot(col2);
2490  T add = 0.0;
2491  for (int k = 0; k<j; ++k) {
2492  if (norms[k] > EPSILON) {
2493  alphat.refCol(k,col2);
2494  add -= G(j,k)*col.dot(col2);
2495  }
2496  }
2497  add += add - G(j,j)*col.nrm2sq();
2498  err2 += add;
2499  }
2500  }
2501  return err2;
2502 }
2503 
2505 template <typename T>
2506 T computeError(const T normX2,
2507  const Matrix<T>& G,const Vector<T>& DtR,const Vector<T>& coeffs,
2508  SpVector<T>& spAlpha) {
2509  coeffs.toSparse(spAlpha);
2510  return normX2 -G.quad(spAlpha)-2*DtR.dot(spAlpha);
2511 };
2512 
2513 /* ******************
2514  * Simultaneous OMP
2515  * *****************/
2516 
2517 template <typename T>
2518 void somp(const Matrix<T>* X, const Matrix<T>& D, SpMatrix<T>* spalpha,
2519  const int Ngroups, const int L, const T eps,const int numThreads) {
2520  somp(X,D,spalpha,Ngroups,L,&eps,false,numThreads);
2521 }
2522 
2523 template <typename T>
2524 void somp(const Matrix<T>* XT, const Matrix<T>& D, SpMatrix<T>* spalphaT,
2525  const int Ngroups, const int LL, const T* eps, const bool adapt,
2526  const int numThreads) {
2527  if (LL <= 0) return;
2528  const int K = D.n();
2529  const int L = MIN(D.m(),MIN(LL,K));
2530 
2531  if (!D.isNormalized()) {
2532  cerr << "Current implementation of OMP does not support non-normalized dictionaries" << endl;
2533  return;
2534  }
2535 
2537  Matrix<T> G;
2538  D.XtX(G);
2539 
2540  int NUM_THREADS=init_omp(numThreads);
2541 
2542  int i;
2543 #pragma omp parallel for private(i)
2544  for (i = 0; i< Ngroups; ++i) {
2545  const Matrix<T>& X = XT[i];
2546  const int M = X.n();
2547  SpMatrix<T>& spalpha = spalphaT[i];
2548  spalpha.clear();
2549  Vector<int> rv;
2550  Matrix<T> vM;
2551  T thrs = adapt ? eps[i] : M*(*eps);
2552  coreSOMP(X,D,G,vM,rv,L,thrs);
2553  spalpha.convert2(vM,rv,K);
2554  }
2555 }
2556 
2557 template <typename T>
2558 void coreSOMP(const Matrix<T>& X, const Matrix<T>& D, const Matrix<T>& G,
2559  Matrix<T>& v,
2560  Vector<int>& r, const int L, const T eps) {
2561  const int K = G.n();
2562  const int n = D.m();
2563  const int M = X.n();
2564 
2565  const bool big_mode = M*K*(n+L) > 2*(M*n*n+K*n*(n+L));
2566  r.resize(L);
2567  r.set(-1);
2568  v.resize(0,X.n());
2569 
2570  if (M == 1) {
2571  Vector<T> scores(K);
2572  Vector<T> norm(K);
2573  Vector<T> tmp(K);
2574  Matrix<T> Un(L,L);
2575  Un.setZeros();
2576  Matrix<T> Undn(K,L);
2577  Matrix<T> Unds(L,L);
2578  Matrix<T> Gs(K,L);
2579  Vector<T> Rdn(K);
2580  Vector<T> Xt(X.rawX(),n);
2581  D.multTrans(Xt,Rdn);
2582  Vector<T> RUn(L);
2583  T normX = Xt.nrm2sq();
2584  T lambda=0;
2585  coreORMP(scores,norm,tmp,Un,Undn,Unds,Gs,Rdn,G,r,RUn,normX,&eps,&L,&lambda);
2586  int count=0;
2587  for (int i = 0; i<L; ++i) {
2588  if (r[i] == -1) break;
2589  ++count;
2590  }
2591  v.resize(count,X.n());
2592  Vector<T> v1(v.rawX(),count);
2593  Vector<T> v2(RUn.rawX(),count);
2594  v1.copy(v2);
2595  return;
2596  }
2597 
2598  Matrix<T> XXtD;
2599  Matrix<T> XtD;
2600  T E;
2601  if (big_mode) {
2602  Matrix<T> XXt;
2603  X.XXt(XXt);
2604  E = XXt.trace();
2605  if (E < eps) return;
2606  XXt.mult(D,XXtD);
2607  } else {
2608  E=X.normFsq();
2609  if (E < eps) return;
2610  X.mult(D,XtD,true);
2611  }
2612 
2613  Matrix<T> A(K,L);
2614  A.setZeros();
2615  Matrix<T> B(L,K);
2616  B.setZeros();
2617  Matrix<T> S(L,L);
2618  S.setZeros();
2619  Matrix<T> Fs(K,L);
2620  Fs.setZeros();
2621  Matrix<T> Gs(K,L);
2622  Gs.setZeros();
2623  Matrix<T> As(L,L);
2624  As.setZeros();
2625 
2626  Vector<T> tmp(K);
2627  Vector<T> e(K);
2628  G.diag(e);
2629  Vector<T> f(K);
2630  if (big_mode) {
2631  for (int i = 0; i<K; ++i) {
2632  Vector<T> di;
2633  D.refCol(i,di);
2634  Vector<T> di2;
2635  XXtD.refCol(i,di2);
2636  f[i]=di.dot(di2);
2637  }
2638  } else {
2639  XtD.norm_2sq_cols(f);
2640  }
2641  Vector<T> c(L);
2642  c.setZeros();
2643  Vector<T> scores(K);
2644 
2646  T* const prAs = As.rawX();
2647  T* const prA = A.rawX();
2648  T* const prS = S.rawX();
2649  T* const prGs = Gs.rawX();
2650  T* const prFs = Fs.rawX();
2651  T* const prB = B.rawX();
2652  T* const pr_c = c.rawX();
2653  T* const pr_tmp = tmp.rawX();
2654 
2655  int j;
2656  for (j = 0; j<L; ++j) {
2657  scores.copy(f);
2658  scores.div(e);
2659  for (int k = 0; k<j; ++k) scores[r[k]]=-1.0;
2660  const int currentInd = scores.max();
2661  const T invNorm=T(1.0)/sqrt(e[currentInd]);
2662  if (invNorm > 1e3) {
2663  j=j-1;
2664  break;
2665  }
2666  r[j]=currentInd;
2667  E -= scores[currentInd];
2668  for (int k = 0; k<j; ++k) prS[j*L+k]=T();
2669  prS[j*L+j]=T(1.0);
2670  for (int k = 0; k<j; ++k) prAs[k*L+j]=prA[k*K+currentInd];
2671 
2673  int iter = invNorm > 1.41 ? 2 : 1;
2674  for (int k = 0; k<iter; ++k) {
2675  for (int l = 0; l<j; ++l) {
2676  T scal = -cblas_dot<T>(j-l+1,prAs+l*L+l,1,prS+j*L+l,1);
2677  cblas_axpy<T>(l+1,scal,prS+l*L,1,prS+j*L,1);
2678  }
2679  }
2680  cblas_scal<T>(j+1,invNorm,prS+j*L,1);
2681 
2682  if (j == L-1 || E <= eps) {
2683  ++j;
2684  break;
2685  }
2686 
2689  Vector<T> Gsj;
2690  Gs.refCol(j,Gsj);
2691  G.copyCol(currentInd,Gsj);
2692  cblas_gemv<T>(CblasColMajor,CblasNoTrans,K,j+1,T(1.0),prGs,K,prS+j*L,1,
2693  T(0.0),prA+j*K,1);
2694  prAs[j*L+j]=prA[j*K+currentInd];
2695  Vector<T> Aj;
2696  A.refCol(j,Aj);
2697  tmp.sqr(Aj);
2698  e.sub(tmp);
2699 
2700  Vector<T> Fsj;
2701  Fs.refCol(j,Fsj);
2702  if (big_mode) {
2703  Vector<T> di;
2704  D.refCol(currentInd,di);
2705  XXtD.multTrans(di,Fsj);
2706  } else {
2707  Vector<T> di;
2708  XtD.refCol(currentInd,di);
2709  XtD.multTrans(di,Fsj);
2710  }
2711  cblas_gemv<T>(CblasColMajor,CblasNoTrans,K,j+1,T(1.0),prFs,K,prS+j*L,1,
2712  T(0.0),prB+j,L);
2713  for (int k = 0; k<j;++k) pr_c[k]=T();
2714  for (int k = 0; k<=j;++k)
2715  cblas_axpy<T>(j,prS[j*L+k],prB+r[k]*L,1,pr_c,1);
2716  f.add(tmp,f[currentInd]*invNorm*invNorm);
2717  if (j > 0) {
2718  cblas_gemv<T>(CblasColMajor,CblasNoTrans,K,j,T(1.0),prA,K,pr_c,1,
2719  T(0.0),pr_tmp,1);
2720  } else {
2721  tmp.setZeros();
2722  }
2723  cblas_axpy<T>(K,T(-1.0),prB+j,L,pr_tmp,1);
2724  tmp.mult(tmp,Aj);
2725  f.add(tmp,T(2.0));
2726  }
2727  A.clear();
2728  B.clear();
2729  Fs.clear();
2730  Gs.clear();
2731  As.clear();
2732 
2733  if (j == 0) return;
2734 
2735  Matrix<T> SSt;
2736  S.upperTriXXt(SSt,j);
2737  Matrix<T> Dg(n,j);
2738  for (int i = 0; i<j;++i) {
2739  Vector<T> Dgi;
2740  Dg.refCol(i,Dgi);
2741  D.copyCol(r[i],Dgi);
2742  }
2743  Matrix<T> SStDt;
2744  SSt.mult(Dg,SStDt,false,true);
2745  SStDt.mult(X,v);
2746 };
2747 
2748 }
2749 
2750 #endif // DECOMP_H
2751 
spams::Matrix::quad
T quad(const Vector< T > &vec1, const SpVector< T > &vec2) const
return vec1&#39;*A*vec2, where vec2 is sparse
Definition: linalg.h:2047
spams::SpMatrix
Sparse Matrix class.
Definition: linalg.h:63
spams::SpMatrix::toFull
void toFull(Matrix< T > &matrix) const
copy the sparse matrix into a dense matrix
Definition: linalg.h:4760
spams
Definition: dag.h:26
spams::Matrix::rank1Update
void rank1Update(const Vector< T > &vec1, const Vector< T > &vec2, const T alpha=1.0)
perform A <- A + alpha*vec1*vec2&#39;
Definition: linalg.h:2384
spams::Vector::sum
T sum() const
returns the sum of the vector
Definition: linalg.h:3344
spams::ist
void ist(const Matrix< T > &X, const Matrix< T > &D, SpMatrix< T > &spalpha, T lambda, constraint_type mode, const int itermax=500, const T tol=0.5, const int numThreads=-1)
Definition: decomp.h:1966
spams::Matrix::norm_2_cols
void norm_2_cols(Vector< T > &norms) const
returns the l2 norms of the columns
Definition: linalg.h:2162
spams::Matrix::meanCol
void meanCol(Vector< T > &mean) const
Compute the mean of the columns.
Definition: linalg.h:2226
CblasLower
Definition: utl_cblas.h:8
CblasRight
Definition: utl_cblas.h:10
spams::computeError
T computeError(const T normX2, const Vector< T > &norms, const Matrix< T > &G, const Matrix< T > &RtD, const Matrix< T > &alphat)
auxiliary function for ist_groupLasso
Definition: decomp.h:2481
utl::vSub
void vSub(int n, T *vecIn, T *vecIn2, T *vecOut)
interface to v*Sub
Definition: utlCoreMKL.h:279
spams::Vector::copy
void copy(const Vector< T > &x)
make a copy of x
Definition: linalg.h:2865
utl::cblas_copy
void cblas_copy(const INTT N, const T *X, const INTT incX, T *Y, const INTT incY)
spams::constraint_type
constraint_type
Definition: decomp.h:88
spams::coreGroupISTConstrained
void coreGroupISTConstrained(const Matrix< T > &G, Matrix< T > &RtD, Matrix< T > &alphat, const T normR, const T eps, const int itermax=500, const T tol=0.5)
Auxiliary function for ist_groupLasso.
Definition: decomp.h:2370
spams::coreISTconstrained
void coreISTconstrained(const AbstractMatrix< T > &G, Vector< T > &DtR, Vector< T > &coeffs, const T normX2, const T thrs, const int itermax=500, const T tol=0.5)
coreIST constrained
Definition: decomp.h:2113
spams::sqr
T sqr(const T x)
template version of the fabs function
spams::Matrix::setn
void setn(const int n)
modify _n
Definition: linalg.h:279
utl::cblas_scal
void cblas_scal(const INTT N, const T alpha, T *X, const INTT incX)
spams::downDateLasso
void downDateLasso(int &j, int &minBasis, T &normX, const bool ols, const bool pos, Vector< T > &Rdn, int *ind, T *coeffs, Vector< T > &sig, Vector< T > &av, Vector< T > &Xdn, Vector< T > &RUn, Matrix< T > &Unm, Matrix< T > &Gsm, Matrix< T > &Gsam, Matrix< T > &Undsm, Matrix< T > &Rm)
Auxiliary functoni for coreLARS (Cholesky downdate)
Definition: decomp.h:1000
spams::Vector::fmaxval
T fmaxval() const
returns the maximum magnitude
Definition: linalg.h:2799
CblasTrans
Definition: utl_cblas.h:6
spams::AbstractMatrix::n
virtual int n() const =0
spams::Vector::add
void add(const Vector< T > &x, const T a=1.0)
A <- A + a*x.
Definition: linalg.h:3029
spams::Matrix::XtX
void XtX(Matrix< T > &XtX) const
XtX = A&#39;*A.
Definition: linalg.h:1959
spams::Matrix::addDiag
void addDiag(const Vector< T > &diag)
Definition: linalg.h:1506
utl::cblas_ger
void cblas_ger(const CBLAS_ORDER order, const INTT M, const INTT N, const T alpha, const T *X, const INTT incX, const T *Y, const INTT incY, T *A, const INTT lda)
utl::E
static constexpr double E
Definition: utlConstants.h:25
spams::Vector::mult
void mult(const Vector< T > &x, const Vector< T > &y)
A <- x .* y.
Definition: linalg.h:3114
spams::Vector::rawX
T * rawX() const
returns a modifiable reference of the data, DANGEROUS
Definition: linalg.h:593
CblasUpper
Definition: utl_cblas.h:8
spams::Matrix::transpose
void transpose(Matrix< T > &trans)
Definition: linalg.h:1489
spams::omp
void omp(const Matrix< T > &X, const Matrix< T > &D, SpMatrix< T > &spalpha, const int *L, const T *eps, const T *lambda, const bool vecL=false, const bool vecEps=false, const bool Lambda=false, const int numThreads=-1, Matrix< T > *path=NULL)
Definition: decomp.h:298
utl::vAdd
void vAdd(int n, T *vecIn, T *vecIn2, T *vecOut)
interface to v*Add
Definition: utlCoreMKL.h:300
spams::Vector::maxval
T maxval() const
returns the maximum value
Definition: linalg.h:2789
spams::Vector::sqr
void sqr(const Vector< T > &x)
A <- x .^ 2.
Definition: linalg.h:3077
spams::coreORMPB
void coreORMPB(Vector< T > &RtD, const AbstractMatrix< T > &G, Vector< int > &ind, Vector< T > &coeffs, T &normX, const int L, const T eps, const T lambda=0)
Auxiliary function of omp.
Definition: decomp.h:498
spams::coreIST
void coreIST(const AbstractMatrix< T > &G, Vector< T > &DtR, Vector< T > &coeffs, const T thrs, const int itermax=500, const T tol=0.5)
coreIST
Definition: decomp.h:2052
spams::Matrix::diag
void diag(Vector< T > &d) const
extract the diagonal
Definition: linalg.h:1985
spams::ProdMatrix::copyCol
void copyCol(const int i, Vector< T > &DtXi) const
compute DtX(:,i)
Definition: linalg.h:5075
spams::lassoReweighted
void lassoReweighted(const Matrix< T > &X, const Matrix< T > &D, SpMatrix< T > &spalpha, int L, const T constraint, constraint_type mode, const bool pos, const T sigma, const int numThreads=-1)
second implementation using matrix inversion lemma
Definition: decomp.h:1113
spams::Vector::setn
void setn(const int n)
Definition: linalg.h:627
spams::coreSOMP
void coreSOMP(const Matrix< T > &X, const Matrix< T > &D, const Matrix< T > &G, Matrix< T > &vM, Vector< int > &rv, const int L, const T eps)
Definition: decomp.h:2558
low
static char low
a few static variables for lapack
Definition: cblas_alt_template.h:20
spams::Vector::alltrue
bool alltrue() const
spams::Vector::copyMask
void copyMask(Vector< T > &out, Vector< bool > &mask) const
extract the rows of a matrix corresponding to a binary mask
Definition: linalg.h:4041
spams::Matrix::copyCol
void copyCol(const int i, Vector< T > &x) const
Copy the column i into x.
Definition: linalg.h:1100
CblasColMajor
Definition: utl_cblas.h:5
spams::Matrix::m
int m() const
Definition: linalg.h:222
CblasNonUnit
Definition: utl_cblas.h:9
spams::somp
void somp(const Matrix< T > *X, const Matrix< T > &D, SpMatrix< T > *spalpha, const int Ngroups, const int L, const T *pr_eps, const bool adapt=false, const int numThreads=-1)
Definition: decomp.h:2524
utils.h
Contains various variables and class timer.
spams::Vector::n
int n() const
returns the size of the vector
Definition: linalg.h:591
spams::Data
Data class, abstract class, useful in the class image.
Definition: linalg.h:130
MIN
#define MIN(a, b)
Definition: utils.h:47
spams::Matrix::refCol
void refCol(int i, Vector< T > &x) const
Reference the column i into the vector x.
Definition: linalg.h:1144
CblasNoTrans
Definition: utl_cblas.h:6
spams::lasso
void lasso(const Matrix< T > &X, const Matrix< T > &D, SpMatrix< T > &spalpha, int L, const T constraint, const T lambda2=0, constraint_type mode=PENALTY, const bool pos=false, const bool ols=false, const int numThreads=-1, Matrix< T > *path=NULL, const int length_path=-1)
Definition: decomp.h:639
spams::Matrix::copyMask
void copyMask(Matrix< T > &out, Vector< bool > &mask) const
extract the rows of a matrix corresponding to a binary mask
Definition: linalg.h:4052
spams::Vector::allfalse
bool allfalse() const
nonUnit
static char nonUnit
Definition: cblas_alt_template.h:22
MAX
#define MAX(a, b)
Definition: utils.h:48
spams::Vector::setZeros
void setZeros()
Set all values to zero.
Definition: linalg.h:2871
spams::Vector::asum
T asum() const
computes the sum of the magnitudes of the vector
Definition: linalg.h:3324
spams::Vector::dot
T dot(const Vector< T > &x) const
returns A&#39;x
Definition: linalg.h:3012
spams::AbstractMatrix::copyCol
virtual void copyCol(const int i, Vector< T > &Xi) const =0
copy X(:,i) into Xi
spams::Matrix::trace
T trace() const
return the trace of the matrix
Definition: linalg.h:2090
spams::abs
T abs(const T x)
template version of the fabs function
spams::Matrix::norm_2sq_cols
void norm_2sq_cols(Vector< T > &norms) const
returns the l2 norms ^2 of the columns
Definition: linalg.h:2204
EPSILON
#define EPSILON
Definition: utils.h:60
spams::Vector
Dense Vector class.
Definition: linalg.h:65
spams::Matrix::isNormalized
bool isNormalized() const
Check wether the columns of the matrix are normalized or not.
Definition: linalg.h:1158
spams::ProdMatrix
Class representing the product of two matrices.
Definition: linalg.h:998
spams::Vector::toSparse
void toSparse(SpVector< T > &vec) const
make a sparse copy
Definition: linalg.h:4025
spams::lasso2
void lasso2(const Matrix< T > &X, const Matrix< T > &D, SpMatrix< T > &spalpha, int L, const T constraint, const T lambda2=0, constraint_type mode=PENALTY, const bool pos=false, const int numThreads=-1, Matrix< T > *path=NULL, const int length_path=-1)
second implementation using matrix inversion lemma
Definition: decomp.h:1468
spams::SpMatrix::clear
void clear()
clear the matrix
Definition: linalg.h:4204
spams::ist_groupLasso
void ist_groupLasso(const Matrix< T > *XT, const Matrix< T > &D, Matrix< T > *alphaT, const int Ngroups, const T lambda, const constraint_type mode, const int itermax=500, const T tol=0.5, const int numThreads=-1)
ist for group Lasso
Definition: decomp.h:2196
spams::Vector::set
void set(const T val)
set each value of the vector to val
Definition: linalg.h:2997
spams::Matrix::fillRow
void fillRow(const Vector< T > &row)
fill the matrix with the row given
Definition: linalg.h:2241
spams::Vector::max
int max() const
returns the index of the largest value
Definition: linalg.h:2761
spams::lassoWeightPreComputed
void lassoWeightPreComputed(const Matrix< T > &X, const Matrix< T > &G, const Matrix< T > &DtR, const Matrix< T > &weights, SpMatrix< T > &spalpha, int L, const T constraint, constraint_type mode, const bool pos, const int numThreads)
Definition: decomp.h:1298
spams::AbstractMatrix
Abstract matrix class.
Definition: linalg.h:188
spams::Matrix::mult
void mult(const Vector< T > &x, Vector< T > &b, const T alpha=1.0, const T beta=0.0) const
perform b = alpha*A*x+beta*b
Definition: linalg.h:1783
spams::Matrix::n
int n() const
Number of columns.
Definition: linalg.h:224
spams::lasso_mask
void lasso_mask(const Matrix< T > &X, const Matrix< T > &D, SpMatrix< T > &spalpha, const Matrix< bool > &mask, int L, const T constraint, const T lambda2=0, constraint_type mode=PENALTY, const bool pos=false, const int numThreads=-1)
second implementation using matrix inversion lemma
Definition: decomp.h:1371
spams::AbstractMatrix::extract_rawCol
virtual void extract_rawCol(const int i, T *Xi) const =0
copy X(:,i) into Xi
SIGN
#define SIGN(a)
Definition: utils.h:49
spams::SpVector
Sparse Vector class.
Definition: linalg.h:67
spams::Vector::nrm2sq
T nrm2sq() const
returns ||A||_2^2
Definition: linalg.h:3007
spams::Matrix::setZeros
void setZeros()
Set all the values to zero.
Definition: linalg.h:1240
spams::Vector::resize
void resize(const int n)
resize the vector
Definition: linalg.h:2876
spams::AbstractMatrix::add_rawCol
virtual void add_rawCol(const int i, T *col, const T a) const =0
compute X(:,i)<- X(:,i)+a*col;
spams::Data::norm_2sq_cols
virtual void norm_2sq_cols(Vector< T > &norms) const
Definition: linalg.h:139
spams::Matrix::XXt
void XXt(Matrix< T > &XXt) const
XXt = A*A&#39;.
Definition: linalg.h:1967
spams::SpMatrix::convert2
void convert2(const Matrix< T > &v, const Vector< int > &r, const int K)
use the data from v, r for _v, _r
Definition: linalg.h:4809
spams::coreLARS2W
void coreLARS2W(Vector< T > &DtR, AbstractMatrix< T > &G, Matrix< T > &Gs, Matrix< T > &Ga, Matrix< T > &invGs, Vector< T > &u, Vector< T > &coeffs, const Vector< T > &weights, Vector< int > &ind, Matrix< T > &work, T &normX, const constraint_type mode, const T constraint, const bool pos=false)
Auxiliary function for lasso.
Definition: decomp.h:1764
spams::Matrix::meanRow
void meanRow(Vector< T > &mean) const
Compute the mean of the rows.
Definition: linalg.h:2233
spams::Matrix::upperTriXXt
void upperTriXXt(Matrix< T > &XXt, const int L) const
XXt = A*A&#39; where A is an upper triangular matrix.
Definition: linalg.h:1975
spams::Vector::fmax
int fmax() const
returns the index of the value with largest magnitude
Definition: linalg.h:2843
spams::SpVector::resize
void resize(const int nzmax)
resizes the vector
Definition: linalg.h:4971
spams::init_omp
static int init_omp(const int numThreads)
Definition: misc.h:264
spams::ProdMatrix::setMatrices
void setMatrices(const Matrix< T > &D, const bool high_memory=true)
set_matrices
Definition: linalg.h:5059
utl::cblas_dot
T cblas_dot(const INTT N, const T *X, const INTT incX, const T *Y, const INTT incY)
spams::coreGroupIST
void coreGroupIST(const Matrix< T > &G, Matrix< T > &RtD, Matrix< T > &alphat, const T thrs, const int itermax=500, const T tol=0.5)
Auxiliary function for ist_groupLasso.
Definition: decomp.h:2290
spams::Matrix::asum
T asum() const
compute the sum of the magnitude of the matrix values
Definition: linalg.h:2085
spams::Matrix::multTrans
void multTrans(const Vector< T > &x, Vector< T > &b, const T alpha=1.0, const T beta=0.0) const
perform b = alpha*A&#39;x + beta*b
Definition: linalg.h:1743
spams::SpMatrix::convert
void convert(const Matrix< T > &v, const Matrix< int > &r, const int K)
use the data from v, r for _v, _r
Definition: linalg.h:4786
spams::Vector::div
void div(const Vector< T > &x)
A <- A ./ x.
Definition: linalg.h:3064
spams::lassoWeight
void lassoWeight(const Matrix< T > &X, const Matrix< T > &D, const Matrix< T > &weights, SpMatrix< T > &spalpha, int L, const T constraint, constraint_type mode, const bool pos, const int numThreads)
Definition: decomp.h:1220
spams::coreLARS2
void coreLARS2(Vector< T > &DtR, const AbstractMatrix< T > &G, Matrix< T > &Gs, Matrix< T > &Ga, Matrix< T > &invGs, Vector< T > &u, Vector< T > &coeffs, Vector< int > &ind, Matrix< T > &work, T &normX, const constraint_type mode, const T constraint, const bool pos=false, T *pr_path=NULL, int length_path=-1)
Auxiliary function for lasso.
Definition: decomp.h:1556
spams::omp_mask
void omp_mask(const Matrix< T > &X, const Matrix< T > &D, SpMatrix< T > &spalpha, const Matrix< bool > &mask, const int *L, const T *eps, const T *lambda, const bool vecL=false, const bool vecEps=false, const bool Lambda=false, const int numThreads=-1, Matrix< T > *path=NULL)
Definition: decomp.h:382
spams::Matrix::resize
void resize(int m, int n)
Resize the matrix.
Definition: linalg.h:1217
spams::Matrix::normFsq
T normFsq() const
return ||A||_F^2
Definition: linalg.h:2110
spams::Matrix::toSparse
void toSparse(SpMatrix< T > &matrix) const
make a sparse copy of the current matrix
Definition: linalg.h:2566
spams::Vector::sub
void sub(const Vector< T > &x)
A <- A - x.
Definition: linalg.h:3052
spams::Matrix
Dense Matrix class.
Definition: linalg.h:61
utl::cblas_nrm2
T cblas_nrm2(const INTT N, const T *X, const INTT incX)
spams::Data::n
virtual int n() const =0
spams::coreLARS
void coreLARS(Vector< T > &Rdn, Vector< T > &Xdn, Vector< T > &A, Vector< T > &u, Vector< T > &sig, Vector< T > &av, Vector< T > &RUn, Matrix< T > &Un, Matrix< T > &Unds, Matrix< T > &Gs, Matrix< T > &Gsa, Matrix< T > &workT, Matrix< T > &R, const AbstractMatrix< T > &G, T &normX, Vector< int > &ind, Vector< T > &coeffs, const T constraint, const bool ols=false, const bool pos=false, constraint_type mode=L1COEFFS, T *path=NULL, int length_path=-1)
Auxiliary function for lasso.
Definition: decomp.h:748
spams::ProdMatrix::addDiag
void addDiag(const T diag)
add something to the diagonal
Definition: linalg.h:5112
spams::coreORMP
void coreORMP(Vector< T > &scores, Vector< T > &norm, Vector< T > &tmp, Matrix< T > &Un, Matrix< T > &Undn, Matrix< T > &Unds, Matrix< T > &Gs, Vector< T > &Rdn, const AbstractMatrix< T > &G, Vector< int > &ind, Vector< T > &RUn, T &normX, const T *eps, const int *L, const T *lambda, T *path=NULL)
Auxiliary function of omp.
Definition: decomp.h:514
spams::Matrix::setm
void setm(const int m)
modify _m
Definition: linalg.h:277
INFINITY
#define INFINITY
Definition: utils.h:62
spams::Matrix::rawX
T * rawX() const
reference a modifiable reference to the data, DANGEROUS
Definition: linalg.h:254
spams::Matrix::clear
void clear()
Clear the matrix.
Definition: linalg.h:1250
Generated by   doxygen 1.8.11