Skip to content

Commit 84f7806

Browse files
dong-haodong-hao
authored
Merge pull request #49 from dong-hao/fix/logic-lineSearchWolfe
Fix line search step logic problem in lineSearchWolfe...
2 parents 225e798 + 1895bca commit 84f7806

2 files changed

Lines changed: 57 additions & 51 deletions

File tree

‎f90/INV/LBFGS.f90‎

Lines changed: 25 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -1122,10 +1122,6 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
11221122
! f'(0) = (df/dm).dot.h
11231123
g_0 = dotProd(grad,h)
11241124

1125-
! setup the lower and upper boundary of alpha
1126-
alpha_l = R_ZERO
1127-
alpha_r = (f_0-(f_0*0.1))/(-g_0)/c2
1128-
11291125
! alpha_1 is the initial step size, which is set in LBFGS
11301126
alpha_1 = alpha
11311127
! compute the trial mHat, f, dHat, eAll, rms
@@ -1179,6 +1175,9 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
11791175
call printf('QUADLS',lambda,alpha,f,mNorm,rms)
11801176
call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile)
11811177
niter = niter + 1
1178+
! use the two evaluated trial points to seed the sectioning bracket
1179+
alpha_l = min(alpha_1,alpha)*0.9 ! give a bit shrinkage to ensure progress
1180+
alpha_r = max(alpha_1,alpha)*1.1 ! give a bit expansion to ensure progress
11821181
! check whether the solution satisfies the sufficient decrease condition
11831182
! Strong Wolfe's condition needs the gradient
11841183
! well, we are going to calculate it anyway - so why don't we do it now?
@@ -1214,14 +1213,15 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
12141213
endif
12151214
endif
12161215
else
1216+
! compute the gradient at the starting guess before using g_1
1217+
call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1)
1218+
g_1 = dotProd(grad, h)
1219+
write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0
1220+
write(*,'(a4,es12.5)') ' g1=',g_1
1221+
write(ioLog,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0
1222+
write(ioLog,'(a4,es12.5)') ' g1=',g_1
12171223
if (f_1 < f_0) then! is the initial making any progress?
12181224
! Test if the initial guess is good for Strong Wolfe condition
1219-
call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1)
1220-
g_1 = dotProd(grad, h)
1221-
write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0
1222-
write(*,'(a4,es12.5)') ' g1=',g_1
1223-
write(ioLog,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0
1224-
write(ioLog,'(a4,es12.5)') ' g1=',g_1
12251225
if ((f_1 <= f_0 + c * alpha_1 * g_0).and.(abs(g_1) <= c2*abs(g_0))) then
12261226
write(*,'(a53)') 'Strong Wolfe Condition satisfied, exiting line search'
12271227
write(ioLog,'(a53)') 'Strong Wolfe Condition satisfied, exiting line search'
@@ -1241,16 +1241,15 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
12411241
call deall_solnVectorMTX(eAll_1)
12421242
return
12431243
endif
1244+
endif
1245+
!ooops, we missed the Strong Wolfe's condition (for one reason or
1246+
!the other
1247+
if ((alpha_r-alpha_l)*g_1<0) then
1248+
! update the left boundary for alpha
1249+
alpha_l = alpha_1
12441250
else
1245-
!ooops, we missed the Strong Wolfe's condition (for one reason or
1246-
!the other
1247-
if ((alpha_r-alpha_l)*g_1<0) then
1248-
! update the left boundary for alpha
1249-
alpha_l = alpha_1
1250-
else
1251-
! update the right boundary for alpha
1252-
alpha_r = alpha_1
1253-
endif
1251+
! update the right boundary for alpha
1252+
alpha_r = alpha_1
12541253
endif
12551254
endif
12561255
! someone has to spread the bad news
@@ -1297,9 +1296,9 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
12971296
if ((b*b-3*a*g_0)<0) then ! failed to fit cubic
12981297
write(*,'(a40)') 'SQRT of negative value in CUBIC INTERP!'
12991298
write(ioLog,'(a40)') 'SQRT of negative value in CUBIC INTERP!'
1300-
write(*,'(a35)') 'using default value to bracket...'
1301-
write(ioLog,'(a35)') 'using default value to bracket...'
1302-
alpha = sqrt(alpha_l*alpha_r)
1299+
write(*,'(a39)') 'using safeguarded midpoint to bracket...'
1300+
write(ioLog,'(a39)') 'using safeguarded midpoint to bracket...'
1301+
alpha = 0.5d0*(alpha_l + alpha_r)
13031302
else
13041303
if (b<=R_ZERO) then ! fit cubic
13051304
alpha = (- b + sqrt(b*b - 3.0*a*g_0))/(3.0*a)
@@ -1311,6 +1310,10 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
13111310
alpha = -g_0/(b+sqrt(b*b - 3.0*a*g_0))
13121311
endif
13131312
endif
1313+
! keep the sectioning step inside the current bracket
1314+
if ((alpha <= alpha_l).or.(alpha >= alpha_r)) then
1315+
alpha = 0.5d0*(alpha_l + alpha_r)
1316+
endif
13141317
! if alpha is too close or too much smaller than its predecessor
13151318
! if ((alpha_j - alpha < eps).or.(alpha < k*alpha_j)) then
13161319
! alpha = alpha_j/TWO ! reset alpha to ensure progress

‎f90/INV/NLCG.f90‎

Lines changed: 32 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -977,7 +977,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
977977
! solution does not satisfy the condition, then backtrack using
978978
! cubic interpolation.
979979
!
980-
! The initial step size is set outside of this routine (in the LBFGS)
980+
! The initial step size is set outside of this routine (in NLCGsolver)
981981
! but these are the good choices (ref. Michael Ferris, Chapter 3, p 59):
982982
! alpha_1 = alpha_{k-1} dotProd(grad_{k-1},h_{k-1})/dotProd(grad_k,h_k})
983983
! or interpolate the quadratic to f(m_{k-1}), f(m_k) and
@@ -1086,11 +1086,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
10861086
! f'(0) = (df/dm).dot.h
10871087
g_0 = dotProd(grad,h)
10881088

1089-
! setup the lower and upper boundary of alpha
1090-
alpha_l = R_ZERO
1091-
alpha_r = (f_0-(f_0*0.1))/(-g_0)/c2
1092-
1093-
! alpha_1 is the initial step size, which is set in LBFGS
1089+
! alpha_1 is the initial step size, which is set in NLCGsolver
10941090
alpha_1 = alpha
10951091
! compute the trial mHat, f, dHat, eAll, rms
10961092
mHat_1 = mHat_0
@@ -1140,6 +1136,9 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
11401136
call printf('QUADLS',lambda,alpha,f,mNorm,rms)
11411137
call printf('QUADLS',lambda,alpha,f,mNorm,rms,logFile)
11421138
niter = niter + 1
1139+
! use the two evaluated trial points to seed the sectioning bracket
1140+
alpha_l = min(alpha_1,alpha)*0.9 ! give a bit shrinkage to ensure progress
1141+
alpha_r = max(alpha_1,alpha)*1.1 ! give a bit expansion to ensure progress
11431142
! check whether the solution satisfies the sufficient decrease condition
11441143
! Strong Wolfe's condition needs the gradient
11451144
! well, we are going to calculate it anyway - so why don't we do it now?
@@ -1172,14 +1171,15 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
11721171
endif
11731172
endif
11741173
else
1174+
! compute the gradient at the starting guess before using g_1
1175+
call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1)
1176+
g_1 = dotProd(grad, h)
1177+
write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0
1178+
write(*,'(a4,es12.5)') ' g1=',g_1
1179+
write(ioLog,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0
1180+
write(ioLog,'(a4,es12.5)') ' g1=',g_1
11751181
if (f_1 < f_0) then! is the initial making any progress?
11761182
! Test if the initial guess is good for Strong Wolfe condition
1177-
call gradient(lambda,d,m0,mHat_1,grad,dHat_1,eAll_1)
1178-
g_1 = dotProd(grad, h)
1179-
write(*,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0
1180-
write(*,'(a4,es12.5)') ' g1=',g_1
1181-
write(ioLog,'(a29,es12.5)',advance='no') ' GRAD: computed, with g0=',g_0
1182-
write(ioLog,'(a4,es12.5)') ' g1=',g_1
11831183
if ((f_1 <= f_0 + c * alpha_1 * g_0).and.(abs(g_1) <= c2*abs(g_0))) then
11841184
write(*,'(a53)') 'Strong Wolfe Condition satisfied, exiting line search'
11851185
write(ioLog,'(a53)') 'Strong Wolfe Condition satisfied, exiting line search'
@@ -1196,21 +1196,20 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
11961196
call deall_solnVectorMTX(eAll_1)
11971197
return
11981198
endif
1199+
endif
1200+
!ooops, we missed the Strong Wolfe's condtion (for one reason or
1201+
!the other
1202+
if ((alpha_r-alpha_l)*g_1<0) then
1203+
! update the left boundary for alpha
1204+
alpha_l = alpha_1
11991205
else
1200-
!ooops, we missed the Strong Wolfe's condtion (for one reason or
1201-
!the other
1202-
if ((alpha_r-alpha_l)*g_1<0) then
1203-
! update the left boundary for alpha
1204-
alpha_l = alpha_1
1205-
else
1206-
! update the right boundary for alpha
1207-
alpha_r = alpha_1
1208-
endif
1206+
! update the right boundary for alpha
1207+
alpha_r = alpha_1
12091208
endif
12101209
endif
12111210
! someone has to spread the bad news
1212-
write(*,'(a50)') '!======Strong Wolfe Condition NOT satisfied!======'
1213-
write(ioLog,'(a50)') '!======Strong Wolfe Condition NOT satisfied!======'
1211+
write(*,'(a47)') '!=====Strong Wolfe Condition NOT satisfied!===='
1212+
write(ioLog,'(a47)') '!=====Strong Wolfe Condition NOT satisfied!===='
12141213

12151214
if (f > f_0) then
12161215
! this should not happen, but in practice it is possible to end up with
@@ -1231,8 +1230,8 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
12311230
write(ioLog,'(a75)') 'Unable to fit a quadratic due to bad gradient estimate, exiting line search'
12321231
else ! /inhales
12331232
! sectioning and fit cubic (initialize)
1234-
write(*,'(a50)') '!========Now sectioning with brute force========'
1235-
write(ioLog,'(a50)') '!========Now sectioning with brute force========'
1233+
write(*,'(a47)') '!=======Now sectioning with brute force========'
1234+
write(ioLog,'(a47)') '!=======Now sectioning with brute force========'
12361235
write(6,*) 'alpha_l=',alpha_l,'alpha_r=',alpha_r
12371236
write(ioLog,*) 'alpha_l=',alpha_l,'alpha_r=',alpha_r
12381237
alpha_i = alpha_1 !START
@@ -1250,16 +1249,20 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
12501249
if ((b*b-3*a*g_0)<0) then ! failed to fit cubic
12511250
write(*,'(a40)') 'SQRT of negative value in CUBIC INTERP!'
12521251
write(ioLog,'(a40)') 'SQRT of negative value in CUBIC INTERP!'
1253-
write(*,'(a35)') 'using default value to bracket...'
1254-
write(ioLog,'(a35)') 'using default value to bracket...'
1255-
alpha = sqrt(alpha_l*alpha_r)
1252+
write(*,'(a39)') 'using safeguarded midpoint to bracket...'
1253+
write(ioLog,'(a39)') 'using safeguarded midpoint to bracket...'
1254+
alpha = 0.5d0*(alpha_l + alpha_r)
12561255
else
12571256
if (b<=R_ZERO) then ! fit cubic
12581257
alpha = (- b + sqrt(b*b - 3.0*a*g_0))/(3.0*a)
12591258
else
12601259
alpha = -g_0/(b+sqrt(b*b - 3.0*a*g_0))
12611260
endif
12621261
endif
1262+
! keep the sectioning step inside the current bracket
1263+
if ((alpha <= alpha_l).or.(alpha >= alpha_r)) then
1264+
alpha = 0.5d0*(alpha_l + alpha_r)
1265+
endif
12631266
! compute the penalty functional
12641267
call linComb(ONE,mHat_0,alpha,h,mHat)
12651268
call func(lambda,d,m0,mHat,f,mNorm,dHat,eAll,rms)
@@ -1299,7 +1302,7 @@ subroutine lineSearchWolfe(lambda,d,m0,h,alpha,mHat,f,grad, &
12991302
f_i = f_j
13001303
alpha_j = alpha
13011304
f_j = f
1302-
if (ibracket >= 5) then
1305+
if (ibracket >= 2) then
13031306
write(*,'(a69)') 'Warning: exiting bracketing since the it has iterated too many times!'
13041307
write(ioLog,'(a69)') 'Warning: exiting bracketing since the it has iterated too many times!'
13051308
if (f < f_1) then

0 commit comments

Comments
 (0)