240323.221234.HKT cobyla: adopt a new way to assess the geometry of i…

…nterpolation set and set `jdrop_geo`, which align better with the other four solvers; consequently `assess_geo` and `setdrop_geo` are removed

fortran/bobyqa/bobyqb.f90

@@ -32,7 +32,7 @@ module bobyqb_mod
 !
 ! Started: February 2022
 !
-! Last Modified: Saturday, March 16, 2024 PM04:44:48
+! Last Modified: Saturday, March 23, 2024 PM10:00:04
 !--------------------------------------------------------------------------------------------------!
 
 implicit none
@@ -324,7 +324,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm
     ! Generate the next trust region step D.
     call trsbox(delta, gopt, hq, pq, sl, su, trtol, xpt(:, kopt), xpt, crvmin, d)
     dnorm = min(delta, norm(d))
-    shortd = (dnorm < HALF * rho)
+    shortd = (dnorm <= HALF * rho)  ! `<=` works better than `<` in case of underflow.
 
     ! Set QRED to the reduction of the quadratic model when the move D is made from XOPT. QRED
     ! should be positive. If it is nonpositive due to rounding errors, we will not take this step.

fortran/cobyla/cobylb.f90

@@ -17,7 +17,7 @@ module cobylb_mod
 !
 ! Started: July 2021
 !
-! Last Modified: Saturday, March 23, 2024 PM05:30:36
+! Last Modified: Saturday, March 23, 2024 PM10:05:37
 !--------------------------------------------------------------------------------------------------!
 
 implicit none
@@ -55,7 +55,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
 use, non_intrinsic :: selectx_mod, only : savefilt, selectx, isbetter
 
 ! Solver-specific modules
-use, non_intrinsic :: geometry_cobyla_mod, only : assess_geo, setdrop_geo, setdrop_tr, geostep
+use, non_intrinsic :: geometry_cobyla_mod, only : setdrop_tr, geostep
 use, non_intrinsic :: initialize_cobyla_mod, only : initxfc, initfilt
 use, non_intrinsic :: trustregion_cobyla_mod, only : trstlp, trrad
 use, non_intrinsic :: update_cobyla_mod, only : updatexfc, updatepole
@@ -343,16 +343,16 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
         exit  ! Better action to take? Geometry step, or simply continue?
     end if
 
-    ! Does the interpolation set have acceptable geometry? It affects IMPROVE_GEO and REDUCE_RHO.
-    adequate_geo = assess_geo(delta, factor_alpha, factor_beta, sim, simi)
+    ! Does the interpolation set have adequate geometry? It affects IMPROVE_GEO and REDUCE_RHO.
+    adequate_geo = all(sum(sim(:, 1:n)**2, dim=1) <= 4.0_RP * delta**2)
 
     ! Calculate the linear approximations to the objective and constraint functions.
     ! N.B.: TRSTLP accesses A mostly by columns, so it is more reasonable to save A instead of A^T.
     ! Zaikun 2023108: According to a test on 2023108, calculating G and A(:, M_LCON+1:M) by solving
     ! the linear systems SIM^T*G = FVAL(1:N)-FVAL(N+1) and SIM^T*A = CONMAT(:, 1:N)-CONMAT(:, N+1)
     ! does not seem to improve or worsen the performance of COBYLA in terms of the number of function
     ! evaluations. The system was solved by SOLVE in LINALG_MOD based on a QR factorization of SIM
-    ! (not necessarily a good algorithm). No preconditioning was used.
+    ! (not necessarily a good algorithm). No preconditioning or scaling was used.
     g = matprod(fval(1:n) - fval(n + 1), simi)
     A(:, 1:m_lcon) = amat
     A(:, m_lcon + 1:m) = transpose(matprod(conmat(m_lcon + 1:m, 1:n) - spread(conmat(m_lcon + 1:m, n + 1), dim=2, ncopies=n), simi))
@@ -368,7 +368,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
     ! TENTH seems to work better than HALF or QUART, especially for linearly constrained problems.
     ! Note that LINCOA has a slightly more sophisticated way of defining SHORTD, taking into account
     ! whether D causes a change to the active set. Should we try the same here?
-    shortd = (dnorm <= TENTH * rho)
+    shortd = (dnorm <= TENTH * rho)  ! `<=` works better than `<` in case of underflow.
 
     ! Predict the change to F (PREREF) and to the constraint violation (PREREC) due to D.
     ! We have the following in precise arithmetic. They may fail to hold due to rounding errors.
@@ -545,26 +545,15 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
     ! we take another geometry step in that case? If no, why should we do it here? Indeed, this
     ! distinction makes no practical difference for CUTEst problems with at most 100 variables
     ! and 5000 constraints, while the algorithm framework is simplified.
-    if (improve_geo .and. .not. assess_geo(delta, factor_alpha, factor_beta, sim, simi)) then
+    if (improve_geo .and. .not. all(sum(sim(:, 1:n)**2, dim=1) <= 4.0_RP * delta**2)) then
         ! Before the geometry step, UPDATEPOLE has been called either implicitly by UPDATEXFC or
         ! explicitly after CPEN is updated, so that SIM(:, N + 1) is the optimal vertex.
 
         ! Decide a vertex to drop from the simplex. It will be replaced with SIM(:, N + 1) + D to
         ! improve acceptability of the simplex. See equations (15) and (16) of the COBYLA paper.
         ! N.B.: COBYLA never sets JDROP_GEO = N + 1.
-        jdrop_geo = setdrop_geo(delta, factor_alpha, factor_beta, sim, simi)
-
-        ! The following JDROP_GEO comes from UOBYQA/NEWUOA/BOBYQA/LINCOA. It performs poorly!
-        !!jdrop_geo = maxloc(sum(sim(:, 1:n)**2, dim=1), dim=1)
-
-        ! JDROP_GEO is between 1 and N unless SIM and SIMI contain NaN, which should not happen
-        ! at this point unless there is a bug. Nevertheless, for robustness, we include the
-        ! following instruction to exit when JDROP_GEO == 0 (if JDROP_GEO does become 0, then
-        ! memory error will occur if we continue, as JDROP_GEO will be used as an index of arrays.)
-        if (jdrop_geo == 0) then
-            info = DAMAGING_ROUNDING
-            exit
-        end if
+        ! The following JDROP_GEO comes from UOBYQA/NEWUOA/BOBYQA/LINCOA.
+        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)

fortran/cobyla/geometry.f90

@@ -8,91 +8,17 @@ module geometry_cobyla_mod
 !
 ! Started: July 2021
 !
-! Last Modified: Tuesday, March 19, 2024 PM08:45:03
+! Last Modified: Saturday, March 23, 2024 PM09:49:28
 !--------------------------------------------------------------------------------------------------!
 
 implicit none
 private
-public :: assess_geo, setdrop_geo, setdrop_tr, geostep
+public :: setdrop_tr, geostep
 
 
 contains
 
 
-function assess_geo(delta, factor_alpha, factor_beta, sim, simi) result(adequate_geo)
-!--------------------------------------------------------------------------------------------------!
-! This function checks if an interpolation set has acceptable geometry as (14) of the COBYLA paper.
-!--------------------------------------------------------------------------------------------------!
-
-! Common modules
-use, non_intrinsic :: consts_mod, only : IK, RP, ONE, TENTH, DEBUGGING
-use, non_intrinsic :: debug_mod, only : assert
-use, non_intrinsic :: infnan_mod, only : is_finite
-use, non_intrinsic :: linalg_mod, only : isinv
-
-implicit none
-
-! Inputs
-real(RP), intent(in) :: delta
-real(RP), intent(in) :: factor_alpha
-real(RP), intent(in) :: factor_beta
-real(RP), intent(in) :: sim(:, :)   ! SIM(N, N+1)
-real(RP), intent(in) :: simi(:, :)  ! SIMI(N, N)
-
-! Outputs
-logical :: adequate_geo
-
-! Local variables
-character(len=*), parameter :: srname = 'ASSESS_GEO'
-integer(IK) :: n
-real(RP) :: veta_scaled(size(sim, 1))
-real(RP) :: vsig_scaled(size(sim, 1))
-real(RP), parameter :: itol = TENTH
-
-! Sizes
-n = int(size(sim, 1), kind(n))
-
-! Preconditions
-if (DEBUGGING) then
-    call assert(n >= 1, 'N >= 1', srname)
-    call assert(size(sim, 1) == n .and. size(sim, 2) == n + 1, 'SIZE(SIM) == [N, N+1]', srname)
-    call assert(all(is_finite(sim)), 'SIM is finite', srname)
-    call assert(all(sum(abs(sim(:, 1:n)), dim=1) > 0), 'SIM(:, 1:N) has no zero column', 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)
-    call assert(isinv(sim(:, 1:n), simi, itol), 'SIMI = SIM(:, 1:N)^{-1}', srname)
-    call assert(delta > 0, 'DELTA > 0', srname)
-    call assert(factor_alpha > 0 .and. factor_alpha < 1, '0 < FACTOR_ALPHA < 1', srname)
-    call assert(factor_beta > 1, 'FACTOR_BETA > 1', srname)
-end if
-
-!====================!
-! Calculation starts !
-!====================!
-
-! Calculate the values of sigma and eta.
-! VETA(J) (1 <= J <= N) is the distance between vertices J and 0 (the best vertex) of the simplex.
-! VSIG(J) is the distance from vertex J to its opposite face of the simplex. Thus VSIG <= VETA.
-! N.B.: What about the distance from vertex N+1 to the its opposite face? Consider the simplex
-! {V_{N+1}, V_{N+1} + L*e_1, ..., V_{N+1} + L*e_N}, where V_{N+1} is vertex N+1, namely the current
-! "best" point, [e_1, ..., e_n] is an orthogonal matrix, and L is a constant in the order of DELTA.
-! This simplex is optimal in the sense that the interpolation system has the minimal condition
-! number, i.e., one. For this simplex, the distance from V_{N+1} to its opposite face is L/SQRT{N}.
-!!vsig = ONE / sqrt(sum(simi**2, dim=2))
-!!veta = sqrt(sum(sim(:, 1:n)**2, dim=1))
-!!adequate_geo = all(vsig >= factor_alpha * delta) .and. all(veta <= factor_beta * delta)
-
-! Zaikun 20240314: We use the scaled versions of VETA and VSIG, as is explained in SETDROP_GEO.
-vsig_scaled = ONE / sqrt(sum((simi * delta)**2, dim=2))
-veta_scaled = sqrt(sum((sim(:, 1:n) / delta)**2, dim=1))
-adequate_geo = all(vsig_scaled >= factor_alpha) .and. all(veta_scaled <= factor_beta)
-
-!====================!
-!  Calculation ends  !
-!====================!
-end function assess_geo
-
-
 function setdrop_tr(ximproved, d, delta, rho, sim, simi) result(jdrop)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine finds (the index) of a current interpolation point to be replaced with the
@@ -268,102 +194,6 @@ function setdrop_tr(ximproved, d, delta, rho, sim, simi) result(jdrop)
 end function setdrop_tr
 
 
-function setdrop_geo(delta, factor_alpha, factor_beta, sim, simi) result(jdrop)
-!--------------------------------------------------------------------------------------------------!
-! This subroutine finds (the index) of a current interpolation point to be replaced with a
-! geometry-improving point. See (15)--(16) of the COBYLA paper.
-! N.B.: COBYLA never sets JDROP = N + 1.
-!--------------------------------------------------------------------------------------------------!
-
-! Common modules
-use, non_intrinsic :: consts_mod, only : IK, RP, ONE, TENTH, DEBUGGING
-use, non_intrinsic :: debug_mod, only : assert
-use, non_intrinsic :: infnan_mod, only : is_nan, is_finite
-use, non_intrinsic :: linalg_mod, only : isinv
-
-implicit none
-
-! Inputs
-real(RP), intent(in) :: delta
-real(RP), intent(in) :: factor_alpha
-real(RP), intent(in) :: factor_beta
-real(RP), intent(in) :: sim(:, :)   ! SIM(N, N+1)
-real(RP), intent(in) :: simi(:, :)  ! SIMI(N, N)
-
-! Outputs
-integer(IK) :: jdrop
-
-! Local variables
-character(len=*), parameter :: srname = 'SETDROP_GEO'
-integer(IK) :: n
-real(RP) :: veta_scaled(size(sim, 1))
-real(RP) :: vsig_scaled(size(sim, 1))
-real(RP), parameter :: itol = TENTH
-
-! Sizes
-n = int(size(sim, 1), kind(n))
-
-! Preconditions
-if (DEBUGGING) then
-    call assert(n >= 1, 'N >= 1', srname)
-    call assert(size(sim, 1) == n .and. size(sim, 2) == n + 1, 'SIZE(SIM) == [N, N+1]', srname)
-    call assert(all(is_finite(sim)), 'SIM is finite', srname)
-    call assert(all(sum(abs(sim(:, 1:n)), dim=1) > 0), 'SIM(:, 1:N) has no zero column', 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)
-    call assert(isinv(sim(:, 1:n), simi, itol), 'SIMI = SIM(:, 1:N)^{-1}', srname)
-    call assert(factor_alpha > 0 .and. factor_alpha < 1, '0 < FACTOR_ALPHA < 1', srname)
-    call assert(factor_beta > 1, 'FACTOR_BETA > 1', srname)
-    call assert(.not. assess_geo(delta, factor_alpha, factor_beta, sim, simi), &
-        & 'The geometry is not adequate when '//srname//' is called', srname)
-end if
-
-!====================!
-! Calculation starts !
-!====================!
-
-! Calculate the values of sigma and eta.
-! VSIG(J) (J=1, .., N) is the Euclidean distance from vertex J to the opposite face of the simplex.
-!!veta = sqrt(sum(sim(:, 1:n)**2, dim=1))
-!!vsig = ONE / sqrt(sum(simi**2, dim=2))
-
-! Zaikun 20240314: We use the scaled versions of VETA and VSIG, which works better if DELTA is small
-! and SIMI contains large values. This prevents an infinite loop that occurred when RP is REAL16
-! (half precision). In that case, VSIG was evaluated to be zero due to rounding, although all entries
-! of VSIG_SCALED ere close to 0.5, and the geometry-improving step turns out to be exactly the same
-! as SIM(:, JDORP), which led to a dead loop.
-veta_scaled = sqrt(sum((sim(:, 1:n) / delta)**2, dim=1))
-vsig_scaled = ONE / sqrt(sum((simi * delta)**2, dim=2))
-
-! Decide which vertex to drop from the simplex. It will be replaced with a new point to improve the
-! acceptability of the simplex. See equations (15) and (16) of the COBYLA paper.
-if (any(veta_scaled > factor_beta)) then
-    jdrop = int(maxloc(veta_scaled, mask=(.not. is_nan(veta_scaled)), dim=1), kind(jdrop))
-    !!MATLAB: [~, jdrop] = max(veta_scaled, [], 'omitnan');
-elseif (any(vsig_scaled < factor_alpha)) then
-    jdrop = int(minloc(vsig_scaled, mask=(.not. is_nan(vsig_scaled)), dim=1), kind(jdrop))
-    !!MATLAB: [~, jdrop] = min(vsig_scaled, [], 'omitnan');
-else
-    ! We arrive here if VSIG_SCALED and VETA_SCALED are all NaN, which can happen due to NaN in SIM
-    ! and SIMI, which should not happen unless there is a bug.
-    jdrop = 0
-end if
-
-! Zaikun 20230202: What if we consider VETA and VSIG together? The following attempts do not work well.
-! !jdrop = maxloc(sum(sim(:, 1:n)**2, dim=1) * sum(simi**2, dim=2)) ! Condition number
-! !jdrop = maxloc(sum(sim(:, 1:n)**2, dim=1)**2 * sum(simi**2, dim=2)) ! Condition number times distance
-
-!====================!
-!  Calculation ends  !
-!====================!
-
-!Postconditions
-if (DEBUGGING) then
-    call assert(jdrop >= 1 .and. jdrop <= n, '1 <= JDROP <= N', srname)
-end if
-end function setdrop_geo
-
-
 function geostep(jdrop, amat, bvec, conmat, cpen, cval, delta, fval, factor_gamma, simi) result(d)
 !--------------------------------------------------------------------------------------------------!
 ! This function calculates a geometry step so that the geometry of the interpolation set is improved

fortran/lincoa/lincob.f90

@@ -15,7 +15,7 @@ module lincob_mod
 !
 ! Started: February 2022
 !
-! Last Modified: Saturday, March 16, 2024 PM04:59:49
+! Last Modified: Saturday, March 23, 2024 PM10:00:56
 !--------------------------------------------------------------------------------------------------!
 
 implicit none
@@ -404,7 +404,8 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b
     ! N.B. The magic number 0.1999 seems to be related to the fact that a linear constraint is
     ! considered nearly active if the point under consideration is within 0.2*DELTA to the boundary
     ! of the constraint. See the subroutine GETACT and Section 3 of Powell (2015) for more details.
-    shortd = ((dnorm < HALF * delta .and. ngetact < 2) .or. dnorm < 0.1999_RP * delta)
+    ! `<=` works better than `<` in case of underflow.
+    shortd = ((dnorm <= HALF * delta .and. ngetact < 2) .or. dnorm <= 0.1999_RP * delta)
     !------------------------------------------------------------------------------------------!
     ! The SHORTD defined above needs NGETACT, which relies on Powell's trust region subproblem
     ! solver. If a different subproblem solver is used, we can take the following SHORTD adopted

fortran/newuoa/newuob.f90

@@ -8,7 +8,7 @@ module newuob_mod
 !
 ! Started: July 2020
 !
-! Last Modified: Saturday, March 16, 2024 PM04:46:01
+! Last Modified: Saturday, March 23, 2024 PM09:59:27
 !--------------------------------------------------------------------------------------------------!
 
 implicit none
@@ -278,7 +278,8 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm
     call trsapp(delta, gopt, hq, pq, trtol, xpt, crvmin, d)
     dnorm = min(delta, norm(d))
     ! SHORTD corresponds to Box 3 of the NEWUOA paper. N.B.: we compare DNORM with RHO, not DELTA.
-    shortd = (dnorm < HALF * rho)  ! HALF seems to work better than TENTH or QUART.
+    ! HALF seems to work better than TENTH or QUART.
+    shortd = (dnorm <= HALF * rho)  ! `<=` works better than `<` in case of underflow.
 
     ! Set QRED to the reduction of the quadratic model when the move D is made from XOPT. QRED
     ! should be positive. If it is nonpositive due to rounding errors, we will not take this step.

fortran/uobyqa/uobyqb.f90

@@ -8,7 +8,7 @@ module uobyqb_mod
 !
 ! Started: February 2022
 !
-! Last Modified: Saturday, March 16, 2024 PM04:45:33
+! Last Modified: Saturday, March 23, 2024 PM10:01:42
 !--------------------------------------------------------------------------------------------------!
 
 implicit none
@@ -274,7 +274,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
 
     ! Check whether D is too short to invoke a function evaluation.
     dnorm = min(delta, norm(d))
-    shortd = (dnorm < HALF * rho)
+    shortd = (dnorm <= HALF * rho)  ! `<=` works better than `<` in case of underflow.
 
     ! Set QRED to the reduction of the quadratic model when the move D is made from XOPT. QRED
     ! should be positive. If it is nonpositive due to rounding errors, we will not take this step.