Skip to content

Commit

Permalink
240324.184105.HKT fortran/cobyla: introduce delbar, remove `factor_…
Browse files Browse the repository at this point in the history 
…alpha`, `factor_beta`, and `factor_gamma`
  • Loading branch information
@zaikunzhang
zaikunzhang committed Mar 24, 2024
1 parent 5ee081e commit 1e7d1f4
Show file tree
Hide file tree
Showing 2 changed files with 21 additions and 42 deletions.
37 changes: 10 additions & 27 deletions fortran/cobyla/cobylb.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ module cobylb_mod
!
! Started: July 2021
!
! Last Modified: Sunday, March 24, 2024 AM02:53:57
! Last Modified: Sunday, March 24, 2024 PM05:17:30
!--------------------------------------------------------------------------------------------------!

implicit none
Expand All @@ -40,8 +40,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et

! Common modules
use, non_intrinsic :: checkexit_mod, only : checkexit
use, non_intrinsic :: consts_mod, only : RP, IK, ZERO, ONE, HALF, QUART, TENTH, EPS, REALMAX, &
& DEBUGGING, MIN_MAXFILT
use, non_intrinsic :: consts_mod, only : RP, IK, ZERO, ONE, HALF, TENTH, EPS, REALMAX, DEBUGGING, MIN_MAXFILT
use, non_intrinsic :: debug_mod, only : assert
use, non_intrinsic :: evaluate_mod, only : evaluate
use, non_intrinsic :: history_mod, only : savehist, rangehist
Expand Down Expand Up @@ -132,6 +131,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
real(RP) :: cpen ! Penalty parameter for constraint in merit function (PARMU in Powell's code)
real(RP) :: cval(size(x) + 1)
real(RP) :: d(size(x))
real(RP) :: delbar
real(RP) :: delta
real(RP) :: distsq(size(x) + 1)
real(RP) :: dnorm
Expand All @@ -158,17 +158,6 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
! 312--314 of Powell's cobylb.f code. Powell's code revises ACTREM to CVAL(N + 1) - CSTRV and PREREM
! to PREREC in this case, which is crucial for feasibility problems.
real(RP), parameter :: cpenmin = EPS
! FACTOR_ALPHA, FACTOR_BETA, and FACTOR_GAMMA are factors that COBYLB uses when managing the
! simplex. Note the following.
! 1. FACTOR_ALPHA < FACTOR_GAMMA < 1 < FACTOR_BETA.
! 2. FACTOR_GAMMA has nothing to do with GAMMA1 and GAMMA2, which are the contracting/expanding
! factors for updating the trust-region radius DELTA.
! 3. Powell used one more factor FACTOR_DELTA = 1.1 (in general, 1 < FACTOR_DELTA <= FACTOR_BETA).
! It had nothing to do with DELTA, which is the trust-region radius. It was used when defining
! JDROP_TR. We use a completely different scheme (see SETDROP_TR), which does not need FACTOR_DELTA.
real(RP), parameter :: factor_alpha = QUART ! The factor alpha in the COBYLA paper
real(RP), parameter :: factor_beta = 2.1_RP ! The factor beta in the COBYLA paper
real(RP), parameter :: factor_gamma = HALF ! The factor gamma in the COBYLA paper

! Sizes
m_lcon = int(size(bvec), kind(m_lcon))
Expand Down Expand Up @@ -202,9 +191,6 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
call assert(size(conhist, 1) == m .and. maxconhist * (maxconhist - maxhist) == 0, &
& 'SIZE(CONHIST, 1) == M, SIZE(CONHIST, 2) == 0 or MAXHIST', srname)
call assert(maxchist * (maxchist - maxhist) == 0, 'SIZE(CHIST) == 0 or MAXHIST', srname)
call assert(factor_alpha > 0 .and. factor_alpha < factor_gamma .and. factor_gamma < 1, &
& '0 < FACTOR_ALPHA < FACTOR_GAMMA < 1', srname)
call assert(factor_beta > 1, 'FACTOR_BETA > 1', srname)
end if

!====================!
Expand Down Expand Up @@ -475,8 +461,8 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et


!----------------------------------------------------------------------------------------------!
! Before the next trust-region iteration, we possibly improve the geometry of simplex or
! reduces RHO according to IMPROVE_GEO and REDUCE_RHO. Now we decide these indicators.
! Before the next trust-region iteration, we possibly improve the geometry of simplex or reduce
! RHO according to IMPROVE_GEO and REDUCE_RHO. Now we decide these indicators.
! N.B.: We must ensure that the algorithm does not set IMPROVE_GEO = TRUE at infinitely many
! consecutive iterations without moving SIM(:, N+1) or reducing RHO. Otherwise, the algorithm
! will get stuck in repetitive invocations of GEOSTEP. This is ensured by the following facts.
Expand Down Expand Up @@ -570,26 +556,23 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
! trial point in the current iteration, which was then replaced with the same geometry-
! improving point in the next iteration, and so on. In this process, the simplex alternated
! between two configurations, neither of which had acceptable geometry. Thus RHO was never
! reduced, leading to infinite cycling. (N.B.: our implementation uses DELTA as the trust
! reduced, leading to infinite cycling. (N.B.: Our implementation uses DELTA as the trust
! region radius, with RHO being its lower bound. When the infinite cycling occurred in this
! test, DELTA = RHO and it could not be reduced due to the requirement that DELTA >= RHO.)
jdrop_geo = int(maxloc(sum(sim(:, 1:n)**2, dim=1), dim=1), kind(jdrop_geo))

! Calculate the geometry step D.
! In NEWUOA, GEOSTEP takes DELBAR = MAX(MIN(TENTH * SQRT(MAXVAL(DISTSQ)), HALF * DELTA), RHO)
! rather than DELTA. This should not be done here, because D should improve the geometry of
! the simplex when SIM(:, JDROP) is replaced with D; the quality of the geometry is defined
! by DELTA instead of DELBAR as in (14) of the COBYLA paper. See GEOSTEP for more detail.
d = geostep(jdrop_geo, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamma, simi)
delbar = HALF * delta
d = geostep(jdrop_geo, amat, bvec, conmat, cpen, cval, delbar, fval, simi)
x = sim(:, n + 1) + d

! If X is close to one of the points in the interpolation set, then we do not evaluate the
! objective and constraints at X, assuming them to have the values at the closest point.
! N.B.:
! 1. If this happens, do NOT include X into the filter, as F and CONSTR are inaccurate.
! 2. In precise arithmetic, the geometry improving step ensures that the distance between X
! and any interpolation point is at least FACTOR_GAMMA*DELTA, yet X may be close to them due
! to rounding. In an experiment with single precision on 20240317, X = SIM(:, N+1) occurred.
! and any interpolation point is at least DELBAR, yet X may be close to them due to
! rounding. In an experiment with single precision on 20240317, X = SIM(:, N+1) occurred.
distsq(n + 1) = sum((x - sim(:, n + 1))**2)
distsq(1:n) = [(sum((x - (sim(:, n + 1) + sim(:, j)))**2), j=1, n)] ! Implied do-loop
!!MATLAB: distsq(1:n) = sum((x - (sim(:,1:n) + sim(:, n+1)))**2, 1) % Implicit expansion
Expand Down
26 changes: 11 additions & 15 deletions fortran/cobyla/geometry.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ module geometry_cobyla_mod
!
! Started: July 2021
!
! Last Modified: Saturday, March 23, 2024 PM09:49:28
! Last Modified: Sunday, March 24, 2024 PM05:13:08
!--------------------------------------------------------------------------------------------------!

implicit none
Expand Down Expand Up @@ -109,7 +109,8 @@ function setdrop_tr(ximproved, d, delta, rho, sim, simi) result(jdrop)
!vsig = ONE / sqrt(sum(simi**2, dim=2))
!sigbar = abs(simid) * vsig
!
!! The following JDROP will overwrite the previous one if its premise holds. FACTOR_DELTA = 1.1.
!! The following JDROP will overwrite the previous one if its premise holds. FACTOR_DELTA = 1.1
!! and FACTOR_ALPHA = 0.25.
!mask = (veta > factor_delta * delta .and. (sigbar >= factor_alpha * delta .or. sigbar >= vsig))
!if (any(mask)) then
! jdrop = int(maxloc(veta, mask=mask, dim=1), kind(jdrop))
Expand Down Expand Up @@ -194,7 +195,7 @@ function setdrop_tr(ximproved, d, delta, rho, sim, simi) result(jdrop)
end function setdrop_tr


function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamma, simi) result(d)
function geostep(jdrop, amat, bvec, conmat, cpen, cval, delbar, fval, simi) result(d)
!--------------------------------------------------------------------------------------------------!
! This function calculates a geometry step so that the geometry of the interpolation set is improved
! when SIM(:, JDRO_GEO) is replaced with SIM(:, N+1) + D. See (15)--(17) of the COBYLA paper.
Expand All @@ -215,8 +216,7 @@ function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamm
real(RP), intent(in) :: conmat(:, :) ! CONMAT(M, N+1)
real(RP), intent(in) :: cpen
real(RP), intent(in) :: cval(:) ! CVAL(N+1)
real(RP), intent(in) :: delta
real(RP), intent(in) :: factor_gamma
real(RP), intent(in) :: delbar
real(RP), intent(in) :: fval(:) ! FVAL(N+1)
real(RP), intent(in) :: simi(:, :) ! SIMI(N, N)

Expand All @@ -242,7 +242,7 @@ function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamm
if (DEBUGGING) then
call assert(m >= m_lcon .and. m >= 0, 'M >= 0', srname)
call assert(n >= 1, 'N >= 1', srname)
call assert(delta > 0, 'DELTA > 0', srname)
call assert(delbar > 0, 'DELBAR > 0', srname)
call assert(cpen > 0, 'CPEN > 0', srname)
call assert(size(simi, 1) == n .and. size(simi, 2) == n, 'SIZE(SIMI) == [N, N]', srname)
call assert(all(is_finite(simi)), 'SIMI is finite', srname)
Expand All @@ -253,20 +253,16 @@ function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamm
call assert(size(cval) == n + 1 .and. .not. any(cval < 0 .or. is_nan(cval) .or. is_posinf(cval)), &
& 'SIZE(CVAL) == NPT and CVAL does not contain negative NaN/+Inf', srname)
call assert(jdrop >= 1 .and. jdrop <= n, '1 <= JDROP <= N', srname)
call assert(factor_gamma > 0 .and. factor_gamma < 1, '0 < FACTOR_GAMMA < 1', srname)
end if

!====================!
! Calculation starts !
!====================!

! SIMI(JDROP, :) is a vector perpendicular to the face of the simplex to the opposite of vertex
! JDROP. Set D to the vector in this direction and with length FACTOR_GAMMA * DELTA. Since
! FACTOR_ALPHA < FACTOR_GAMMA < FACTOR_BETA, D improves the geometry of the simplex as per (14) of
! the COBYLA paper. This also explains why this subroutine does not replace DELTA with
! DELBAR = MAX(MIN(TENTH * SQRT(MAXVAL(DISTSQ)), HALF * DELTA), RHO) as in NEWUOA.
! JDROP. Set D to the vector in this direction and with length DELBAR.
d = simi(jdrop, :)
d = factor_gamma * delta * (d / norm(d))
d = delbar * (d / norm(d))

! The code below chooses the direction of D according to an approximation of the merit function.
! See (17) of the COBYLA paper and line 225 of Powell's cobylb.f.
Expand All @@ -293,10 +289,10 @@ function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamm
! Postconditions
if (DEBUGGING) then
call assert(size(d) == n .and. all(is_finite(d)), 'SIZE(D) == N, D is finite', srname)
! In theory, ||S|| == FACTOR_GAMMA*DELTA, which may be false due to rounding, but not too far.
! In theory, ||S|| == DELBAR, which may be false due to rounding, but not too far.
! It is crucial to ensure that the geometry step is nonzero, which holds in theory.
call assert(norm(d) > 0.9_RP * factor_gamma * delta .and. norm(d) <= 1.1_RP * factor_gamma * delta, &
& '||D|| == FACTOR_GAMMA*DELTA', srname)
call assert(norm(d) > 0.9_RP * delbar .and. norm(d) <= 1.1_RP * delbar, &
& '||D|| == DELBAR', srname)
end if
end function geostep

Expand Down

0 comments on commit 1e7d1f4

Please sign in to comment.