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

removal of redundant memory copies, stopping when grad norm too low

This commit is contained in:
Sarod Yatawatta 2018-06-12 12:47:06 +02:00
parent 67f71e0ac5
commit a415926866
3 changed files with 33 additions and 17 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);

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

@ -460,8 +460,8 @@ cubic_interp(
} else {
/* evaluate function for this root */
//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+(z0)*pk */
my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(z0)*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);

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

@ -314,8 +314,8 @@ cubic_interp(
} else {
/* evaluate function for this root */
//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+(z0)*pk */
my_daxpy(m,pk,-b+step+a+z0*(b-a),xp); /* xp<=xp+(z0)*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);