p61743098
/
sagecal
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Merge pull request #49 from nlesc-dirac/dev 

streamlining LBFGS routines
This commit is contained in:
Hanno Spreeuw 2018-06-13 10:10:13 +02:00 committed by GitHub
parent b38add6d20 a415926866
commit b72644b972
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 79 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/lib/Dirac/lbfgs.c
View File

@ -351,8 +351,9 @@ cubic_interp(
// f0d=(p01*p01-p02*p02)/(2.0*step);
f0d=(p01-p02)/(2.0*step);
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,b,xp); /* xp<=xp+(b)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(a-step)*pk */
//my_daxpy(m,pk,b,xp); /* xp<=xp+(b)*pk */
my_daxpy(m,pk,-a+step+b,xp); /* xp<=xp+(b)*pk */
// func(xp,x,m,n,adata);
// my_daxpy(n,xo,-1.0,x);
// f1=my_dnrm2(n,x);
@ -416,8 +417,9 @@ cubic_interp(
fz0=f0+f1;
} else {
/* evaluate function for this root */
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(z0)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(b-step)*pk */
//my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
// func(xp,x,m,n,adata);
// my_daxpy(n,xo,-1.0,x);
// fz0=my_dnrm2(n,x);
@ -510,8 +512,9 @@ linesearch_zoom(
/* evaluate phi(aj) */
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,aj,xp); /* xp<=xp+(alphaj)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp = xk+alphaj*pk */
//my_daxpy(m,pk,aj,xp); /* xp<=xp+(aj)*pk */
my_daxpy(m,pk,-alphaj+aj,xp); /* xp<=xp+(aj)*pk */
// func(xp,x,m,n,adata);
/* calculate x<=x-xo */
// my_daxpy(n,xo,-1.0,x);
@ -530,8 +533,9 @@ linesearch_zoom(
bj=alphaj; /* aj unchanged */
} else {
/* evaluate grad(alphaj) */
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp = xk+aj*pk */
//my_daxpy(m,pk,alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
my_daxpy(m,pk,-aj+alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
// func(xp,x,m,n,adata);
/* calculate x<=x-xo */
// my_daxpy(n,xo,-1.0,x);
@ -696,6 +700,16 @@ linesearch(
#ifdef DEBUG
printf("cost=%lf grad=%lf mu=%lf, alpha1=%lf\n",phi_0,gphi_0,mu,alpha1);
#endif
/* catch if not finite (deltaphi=0 or nan) */
if (!isnormal(mu)) {
free(x);
free(xp);
#ifdef DEBUG
printf("line interval too small\n");
#endif
return mu;
}
ci=1;
alphai=alpha1; /* initial value for alpha(i) : check if 0<alphai<=mu */
@ -737,8 +751,9 @@ linesearch(
}
/* evaluate grad(phi(alpha(i))) */
my_dcopy(m,xk,1,xp,1); /* NOT NEEDED here?? xp<=xk */
my_daxpy(m,pk,alphai+step,xp); /* xp<=xp+(alphai+step)*pk */
//my_dcopy(m,xk,1,xp,1); /* NOT NEEDED here because already xp=xk+alphai*pk */
//my_daxpy(m,pk,alphai+step,xp); /* xp<=xp+(alphai+step)*pk */
my_daxpy(m,pk,step,xp); /* xp<=xp+(alphai+step)*pk */
// func(xp,x,m,n,adata);
/* calculate x<=x-xo */
// my_daxpy(n,xo,-1.0,x);
@ -998,7 +1013,7 @@ lbfgs_fit_common(
cm=0;
ci=0;
while (ck<itmax) {
while (ck<itmax && isnormal(gradnrm) && gradnrm>CLM_STOP_THRESH) {
/* mult with hessian pk=-H_k*gk */
if (ck<M) {
mult_hessian(m,pk,gk,s,y,rho,ck,ci);
@ -1011,10 +1026,11 @@ lbfgs_fit_common(
/* parameters alpha1=10.0,sigma=0.1, rho=0.01, t1=9, t2=0.1, t3=0.5 */
/* FIXME: update paramters for GPU gradient */
alphak=linesearch(func,xk,pk,10.0,0.1,0.01,9,0.1,0.5,x,m,n,step,adata,&tp, &tpg);
/* check if step size is too small, then stop */
if (fabs(alphak)<CLM_STOP_THRESH) {
/* check if step size is too small, or nan, then stop */
if (!isnormal(alphak) || fabs(alphak)<CLM_EPSILON) {
break;
}
/* update parameters xk1=xk+alpha_k *pk */
my_dcopy(m,xk,1,xk1,1);
my_daxpy(m,pk,alphak,xk1);

31
src/lib/Dirac/lbfgs_nocuda.c
View File

@ -421,8 +421,9 @@ cubic_interp(
p02=my_dnrm2(n,x);
f0d=(p01*p01-p02*p02)/(2.0*step);
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,b,xp); /* xp<=xp+(b)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(a-step)*pk */
//my_daxpy(m,pk,b,xp); /* xp<=xp+(b)*pk */
my_daxpy(m,pk,-a+step+b,xp); /* xp<=xp+(b)*pk */
func(xp,x,m,n,adata);
my_daxpy(n,xo,-1.0,x);
f1=my_dnrm2(n,x);
@ -458,8 +459,9 @@ cubic_interp(
fz0=f0+f1;
} else {
/* evaluate function for this root */
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(z0)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(b-step)*pk */
//my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
func(xp,x,m,n,adata);
my_daxpy(n,xo,-1.0,x);
fz0=my_dnrm2(n,x);
@ -540,8 +542,9 @@ linesearch_zoom(
phi_j*=phi_j;
/* evaluate phi(aj) */
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,aj,xp); /* xp<=xp+(alphaj)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp = xk+alphaj*pk */
//my_daxpy(m,pk,aj,xp); /* xp<=xp+(aj)*pk */
my_daxpy(m,pk,-alphaj+aj,xp); /* xp<=xp+(aj)*pk */
func(xp,x,m,n,adata);
/* calculate x<=x-xo */
my_daxpy(n,xo,-1.0,x);
@ -553,8 +556,9 @@ linesearch_zoom(
bj=alphaj; /* aj unchanged */
} else {
/* evaluate grad(alphaj) */
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp = xk+aj*pk */
//my_daxpy(m,pk,alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
my_daxpy(m,pk,-aj+alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
func(xp,x,m,n,adata);
/* calculate x<=x-xo */
my_daxpy(n,xo,-1.0,x);
@ -679,6 +683,7 @@ linesearch(
#endif
/* catch if not finite (deltaphi=0 or nan) */
if (!isnormal(mu)) {
free(x);
free(xp);
#ifdef DEBUG
printf("line interval too small\n");
@ -718,8 +723,9 @@ linesearch(
}
/* evaluate grad(phi(alpha(i))) */
my_dcopy(m,xk,1,xp,1); /* NOT NEEDED here?? xp<=xk */
my_daxpy(m,pk,alphai+step,xp); /* xp<=xp+(alphai+step)*pk */
//my_dcopy(m,xk,1,xp,1); /* NOT NEEDED here because already xp=xk+alphai*pk */
//my_daxpy(m,pk,alphai+step,xp); /* xp<=xp+(alphai+step)*pk */
my_daxpy(m,pk,step,xp); /* xp<=xp+(alphai+step)*pk */
func(xp,x,m,n,adata);
/* calculate x<=x-xo */
my_daxpy(n,xo,-1.0,x);
@ -875,10 +881,11 @@ lbfgs_fit(
/* parameters c1=1e-4 c2=0.9, alpha1=1.0, alphamax=10.0, step (for alpha)=1e-4*/
//alphak=linesearch_nw(func,xk,pk,1.0,10.0,1e-4,0.9,x,m,n,1e-4,adata);
//alphak=1.0;
/* check if step size is too small, then stop */
if (fabs(alphak)<CLM_EPSILON) {
/* check if step size is too small, or nan, then stop */
if (!isnormal(alphak) || fabs(alphak)<CLM_EPSILON) {
break;
}
/* update parameters xk1=xk+alpha_k *pk */
my_dcopy(m,xk,1,xk1,1);
my_daxpy(m,pk,alphak,xk1);

42
src/lib/Dirac/robust_lbfgs_nocuda.c
View File

@ -272,8 +272,9 @@ cubic_interp(
p02=func_robust(x,xo,n,dp->robust_nu,dp->Nt);
f0d=(p01-p02)/(2.0*step);
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,b,xp); /* xp<=xp+(b)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(a-step)*pk */
//my_daxpy(m,pk,b,xp); /* xp<=xp+(b)*pk */
my_daxpy(m,pk,-a+step+b,xp); /* xp<=xp+(b)*pk */
func(xp,x,m,n,adata);
//my_daxpy(n,xo,-1.0,x);
//f1=my_dnrm2(n,x);
@ -312,8 +313,9 @@ cubic_interp(
fz0=f0+f1;
} else {
/* evaluate function for this root */
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(z0)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp=xk+(b-step)*pk */
//my_daxpy(m,pk,a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(a+z0(b-a))*pk */
func(xp,x,m,n,adata);
//my_daxpy(n,xo,-1.0,x);
//fz0=my_dnrm2(n,x);
@ -394,8 +396,9 @@ linesearch_zoom(
phi_j=func_robust(x,xo,n,dp->robust_nu,dp->Nt);
/* evaluate phi(aj) */
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,aj,xp); /* xp<=xp+(alphaj)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp = xk+alphaj*pk */
//my_daxpy(m,pk,aj,xp); /* xp<=xp+(aj)*pk */
my_daxpy(m,pk,-alphaj+aj,xp); /* xp<=xp+(aj)*pk */
func(xp,x,m,n,adata);
/* calculate x<=x-xo */
//my_daxpy(n,xo,-1.0,x);
@ -410,8 +413,9 @@ linesearch_zoom(
bj=alphaj; /* aj unchanged */
} else {
/* evaluate grad(alphaj) */
my_dcopy(m,xk,1,xp,1); /* xp<=xk */
my_daxpy(m,pk,alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
//my_dcopy(m,xk,1,xp,1); /* xp<=xk: NOT NEEDED because already xp = xk+aj*pk */
//my_daxpy(m,pk,alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
my_daxpy(m,pk,-aj+alphaj+step,xp); /* xp<=xp+(alphaj+step)*pk */
func(xp,x,m,n,adata);
/* calculate x<=x-xo */
//my_daxpy(n,xo,-1.0,x);
@ -546,6 +550,16 @@ linesearch(
#ifdef DEBUG
printf("cost=%lf grad=%lf mu=%lf, alpha1=%lf\n",phi_0,gphi_0,mu,alpha1);
#endif
/* catch if not finite (deltaphi=0 or nan) */
if (!isnormal(mu)) {
free(x);
free(xp);
#ifdef DEBUG
printf("line interval too small\n");
#endif
return mu;
}
ci=1;
alphai=alpha1; /* initial value for alpha(i) : check if 0<alphai<=mu */
@ -580,8 +594,9 @@ linesearch(
}
/* evaluate grad(phi(alpha(i))) */
my_dcopy(m,xk,1,xp,1); /* NOT NEEDED here?? xp<=xk */
my_daxpy(m,pk,alphai+step,xp); /* xp<=xp+(alphai+step)*pk */
//my_dcopy(m,xk,1,xp,1); /* NOT NEEDED here because already xp=xk+alphai*pk */
//my_daxpy(m,pk,alphai+step,xp); /* xp<=xp+(alphai+step)*pk */
my_daxpy(m,pk,step,xp); /* xp<=xp+(alphai+step)*pk */
func(xp,x,m,n,adata);
/* calculate x<=x-xo */
//my_daxpy(n,xo,-1.0,x);
@ -993,7 +1008,7 @@ lbfgs_fit_robust(
cm=0;
ci=0;
while (ck<itmax) {
while (ck<itmax && isnormal(gradnrm) && gradnrm>CLM_STOP_THRESH) {
/* mult with hessian pk=-H_k*gk */
if (ck<M) {
mult_hessian(m,pk,gk,s,y,rho,ck,ci);
@ -1008,6 +1023,11 @@ lbfgs_fit_robust(
/* parameters c1=1e-4 c2=0.9, alpha1=1.0, alphamax=10.0, step (for alpha)=1e-4*/
//alphak=linesearch_nw(func_robust,xk,pk,1.0,10.0,1e-4,0.9,x,m,n,1e-4,adata);
//alphak=1.0;
/* check if step size is too small, or nan, then stop */
if (!isnormal(alphak) || fabs(alphak)<CLM_EPSILON) {
break;
}
/* update parameters xk1=xk+alpha_k *pk */
my_dcopy(m,xk,1,xk1,1);
my_daxpy(m,pk,alphak,xk1);