This subroutine seeks the least value of a function of many variables, by applying a trust region method that forms quadratic models by interpolation. There is usually some freedom in the interpolation conditions, which is taken up by minimizing the Frobenius norm of the change to the second derivative of the model, beginning with the zero matrix. The values of the variables are constrained by upper and lower bounds.
In addition to providing CALFUN, an initial vector of variables and
the lower and upper bounds, the user has to set the values of the parameters
RHOBEG, RHOEND and NPT. After scaling the individual variables
if necessary, so that the magnitudes of their expected changes are similar,
RHOBEG is the initial steplength for changes to the variables, a reasonable choice
being the mesh size of a coarse grid search. Further, RHOEND should be suitable for
a search on a very fine grid. Typically, the software calculates a vector
of variables that is within distance 10*RHOEND of a local minimum. Another
consideration is that every trial vector of variables is forced to satisfy
the lower and upper bounds, but there has to be room to make a search in all
directions. Therefore an error return occurs if the difference between the
bounds on any variable is less than 2*RHOBEG. The parameter NPT specifies
the number of interpolation conditions on each quadratic model, the value
NPT=2*N+1 being recommended for a start, where N is the number of
variables. It is often worthwhile to try other choices too, but much larger values
tend to be inefficient, because the amount of routine work of each iteration is
of magnitude NPT**2, and because the achievement of adequate accuracy in some
matrix calculations becomes more difficult. Some excellent numerical results
have been found in the case NPT=N+6 even with more than 100 variables.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n |
number of variables (must be at least two) |
||
| integer, | intent(in) | :: | npt |
number of interpolation conditions. Its value must be in the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not recommended. |
||
| real(kind=wp), | intent(inout), | dimension(:) | :: | x |
Initial values of the variables must be set in X(1),X(2),...,X(N). They will be changed to the values that give the least calculated F. |
|
| real(kind=wp), | intent(in), | dimension(:) | :: | xl |
lower bounds on x. The construction of quadratic models requires XL(I) to be strictly less than XU(I) for each I. Further, the contribution to a model from changes to the I-th variable is damaged severely by rounding errors if XU(I)-XL(I) is too small. |
|
| real(kind=wp), | intent(in), | dimension(:) | :: | xu |
upper bounds on x. The construction of quadratic models requires XL(I) to be strictly less than XU(I) for each I. Further, the contribution to a model from changes to the I-th variable is damaged severely by rounding errors if XU(I)-XL(I) is too small. |
|
| real(kind=wp), | intent(in) | :: | rhobeg |
RHOBEG must be set to the initial value of a trust region radius. It must be positive, and typically should be about one tenth of the greatest expected change to a variable. An error return occurs if any of the differences XU(I)-XL(I), I=1,...,N, is less than 2*RHOBEG. |
||
| real(kind=wp), | intent(in) | :: | rhoend |
RHOEND must be set to the final value of a trust region radius. It must be positive with RHOEND no greater than RHOBEG. Typically, RHOEND should indicate the accuracy that is required in the final values of the variables. |
||
| integer, | intent(in) | :: | iprint |
IPRINT should be set to 0, 1, 2 or 3, which controls the amount of printing. Specifically, there is no output if IPRINT=0 and there is output only at the return if IPRINT=1. Otherwise, each new value of RHO is printed, with the best vector of variables so far and the corresponding value of the objective function. Further, each new value of F with its variables are output if IPRINT=3. |
||
| integer, | intent(in) | :: | maxfun |
an upper bound on the number of calls of CALFUN. |
||
| procedure(func) | :: | calfun |
SUBROUTINE CALFUN (N,X,F) has to be provided by the user. It must set F to the value of the objective function for the current values of the variables X(1),X(2),...,X(N), which are generated automatically in a way that satisfies the bounds given in XL and XU. |
subroutine bobyqa (n, npt, x, xl, xu, rhobeg, rhoend, iprint, maxfun, calfun) implicit none integer,intent(in) :: n !! number of variables (must be at least two) integer,intent(in) :: npt !! number of interpolation conditions. Its value must be in !! the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not !! recommended. real(wp),dimension(:),intent(inout) :: x !! Initial values of the variables must be set in X(1),X(2),...,X(N). They !! will be changed to the values that give the least calculated F. real(wp),dimension(:),intent(in) :: xl !! lower bounds on x. The construction of quadratic models !! requires XL(I) to be strictly less than XU(I) for each I. Further, !! the contribution to a model from changes to the I-th variable is !! damaged severely by rounding errors if XU(I)-XL(I) is too small. real(wp),dimension(:),intent(in) :: xu !! upper bounds on x. The construction of quadratic models !! requires XL(I) to be strictly less than XU(I) for each I. Further, !! the contribution to a model from changes to the I-th variable is !! damaged severely by rounding errors if XU(I)-XL(I) is too small. real(wp),intent(in) :: rhobeg !! RHOBEG must be set to the initial value of a trust region radius. !! It must be positive, and typically should be about one tenth of the greatest !! expected change to a variable. An error return occurs if any of !! the differences XU(I)-XL(I), I=1,...,N, is less than 2*RHOBEG. real(wp),intent(in) :: rhoend !! RHOEND must be set to the final value of a trust !! region radius. It must be positive with RHOEND no greater than !! RHOBEG. Typically, RHOEND should indicate the !! accuracy that is required in the final values of the variables. integer,intent(in) :: iprint !! IPRINT should be set to 0, 1, 2 or 3, which controls the !! amount of printing. Specifically, there is no output if IPRINT=0 and !! there is output only at the return if IPRINT=1. Otherwise, each new !! value of RHO is printed, with the best vector of variables so far and !! the corresponding value of the objective function. Further, each new !! value of F with its variables are output if IPRINT=3. integer,intent(in) :: maxfun !! an upper bound on the number of calls of CALFUN. procedure (func) :: calfun !! SUBROUTINE CALFUN (N,X,F) has to be provided by the user. It must set !! F to the value of the objective function for the current values of the !! variables X(1),X(2),...,X(N), which are generated automatically in a !! way that satisfies the bounds given in XL and XU. integer :: ibmat,id,ifv,igo,ihq,ipq,isl,isu,ivl,iw,ixa,& ixb,ixn,ixo,ixp,izmat,j,jsl,jsu,ndim,np real(wp),dimension(:),allocatable :: w real(wp) :: temp real(wp),parameter :: zero = 0.0_wp ! The array W will be used for working space. allocate( w((NPT+5)*(NPT+N)+3*N*(N+5)/2) ) ! ! Return if the value of NPT is unacceptable. ! np = n + 1 if (npt < n+2 .or. npt > ((n+2)*np)/2) then write(*,'(/4X,A)') & 'Return from BOBYQA because NPT is not in the required interval' return end if ! ! Partition the working space array, so that different parts of it can ! be treated separately during the calculation of BOBYQB. The partition ! requires the first (NPT+2)*(NPT+N)+3*N*(N+5)/2 elements of W plus the ! space that is taken by the last array in the argument list of BOBYQB. ! ndim = npt + n ixb = 1 ixp = ixb + n ifv = ixp + n * npt ixo = ifv + npt igo = ixo + n ihq = igo + n ipq = ihq + (n*np) / 2 ibmat = ipq + npt izmat = ibmat + ndim * n isl = izmat + npt * (npt-np) isu = isl + n ixn = isu + n ixa = ixn + n id = ixa + n ivl = id + n iw = ivl + ndim ! ! Return if there is insufficient space between the bounds. Modify the ! initial X if necessary in order to avoid conflicts between the bounds ! and the construction of the first quadratic model. The lower and upper ! bounds on moves from the updated X are set now, in the ISL and ISU ! partitions of W, in order to provide useful and exact information about ! components of X that become within distance RHOBEG from their bounds. ! do j = 1, n temp = xu (j) - xl (j) if (temp < rhobeg+rhobeg) then write(*,'(/4X,A)') & 'Return from BOBYQA because one of the differences '//& 'XU(I)-XL(I) is less than 2*RHOBEG.' return end if jsl = isl + j - 1 jsu = jsl + n w (jsl) = xl (j) - x (j) w (jsu) = xu (j) - x (j) if (w(jsl) >=-rhobeg) then if (w(jsl) >= zero) then x (j) = xl (j) w (jsl) = zero w (jsu) = temp else x (j) = xl (j) + rhobeg w (jsl) = - rhobeg w (jsu) = max (xu(j)-x(j), rhobeg) end if else if (w(jsu) <= rhobeg) then if (w(jsu) <= zero) then x (j) = xu (j) w (jsl) = - temp w (jsu) = zero else x (j) = xu (j) - rhobeg w (jsl) = min (xl(j)-x(j),-rhobeg) w (jsu) = rhobeg end if end if end do ! ! Make the call of BOBYQB. ! call bobyqb (n, npt, x, xl, xu, rhobeg, rhoend, iprint, maxfun, w(ixb), w(ixp), & w(ifv), w(ixo), w(igo), w(ihq), w(ipq), w(ibmat), w(izmat), ndim, w(isl), & w(isu), w(ixn), w(ixa), w(id), w(ivl), w(iw), calfun) deallocate(w) end subroutine bobyqa