BCODH

Minimizes a function of N variables subject to bounds on the variables using a modified Newton method and a finite-difference Hessian.

Required Arguments

FCN — User-supplied subroutine to evaluate the function to be minimized. The usage is

CALL FCN (N, X, F), where

CALL FCN (N, X, F), where

N – Length of X. (Input)

X – Vector of length N at which point the function is evaluated. (Input)

X should not be changed by FCN.

X should not be changed by FCN.

F – The computed function value at the point X. (Output)

FCN must be declared EXTERNAL in the calling program.

GRAD — User-supplied subroutine to compute the gradient at the point X. The usage is CALL GRAD (N, X, G), where

N – Length of X and G. (Input)

X – Vector of length N at which point the gradient is evaluated. (Input)

X should not be changed by GRAD.

X should not be changed by GRAD.

G – The gradient evaluated at the point X. (Output)

GRAD must be declared EXTERNAL in the calling program.

IBTYPE — Scalar indicating the types of bounds on variables. (Input)

IBTYPE Action

0 User will supply all the bounds.

1 All variables are nonnegative.

2 All variables are nonpositive.

3 User supplies only the bounds on 1st variable, all other variables will have the same bounds.

XLB — Vector of length N containing the lower bounds on the variables. (Input)

XUB — Vector of length N containing the upper bounds on the variables. (Input)

X — Vector of length N containing the computed solution. (Output)

Optional Arguments

N — Dimension of the problem. (Input)

Default: N = SIZE (X,1).

Default: N = SIZE (X,1).

XGUESS — Vector of length N containing the initial guess of the minimum. (Input)

Default: XGUESS = 0.0.

Default: XGUESS = 0.0.

XSCALE — Vector of length N containing the diagonal scaling matrix for the variables. (Input)

XSCALE is used mainly in scaling the gradient and the distance between two points. In the absence of other information, set all entries to 1.0.

Default: XSCALE = 1.0.

XSCALE is used mainly in scaling the gradient and the distance between two points. In the absence of other information, set all entries to 1.0.

Default: XSCALE = 1.0.

FSCALE — Scalar containing the function scaling. (Input)

FSCALE is used mainly in scaling the gradient. In the absence of other information, set FSCALE to 1.0.

Default: FSCALE = 1.0.

FSCALE is used mainly in scaling the gradient. In the absence of other information, set FSCALE to 1.0.

Default: FSCALE = 1.0.

IPARAM — Parameter vector of length 7. (Input/Output)

Set IPARAM(1) to zero for default values of IPARAM and RPARAM. See Comment 4.

Default: IPARAM = 0.

Set IPARAM(1) to zero for default values of IPARAM and RPARAM. See Comment 4.

Default: IPARAM = 0.

RPARAM — Parameter vector of length 7. (Input/Output)

See Comment 4.

See Comment 4.

FVALUE — Scalar containing the value of the function at the computed solution. (Output)

FORTRAN 90 Interface

Generic: CALL BCODH (FCN, GRAD, IBTYPE, XLB, XUB, X [, …])

Specific: The specific interface names are S_BCODH and D_BCODH.

FORTRAN 77 Interface

Single: CALL BCODH (FCN, GRAD, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FVALUE)

Double: The double precision name is DBCODH.

Description

The routine BCODH uses a modified Newton method and an active set strategy to solve minimization problems subject to simple bounds on the variables. The problem is stated as

subject to l ≤ x ≤ u

From a given starting point xc, an active set IA, which contains the indices of the variables at their bounds, is built. A variable is called a “free variable” if it is not in the active set. The routine then computes the search direction for the free variables according to the formula

d = ‑H-1 gc

where H is the Hessian and gc is the gradient evaluated at xc; both are computed with respect to the free variables. The search direction for the variables in IA is set to zero. A line search is used to find a new point xn ,

xn = xc + λd, λ ∈ (0, 1]

such that

f (xn) ≤ f (xc) + αgT d, α∈ (0, 0.5)

Finally, the optimality conditions

∥g(xi)∥ ≤ ɛ, li < xi< ui

g(xi) < 0, xi = ui

g(xi) > 0, xi = li

are checked where ɛ is a gradient tolerance. When optimality is not achieved, another search direction is computed to begin the next iteration. This process is repeated until the optimality criterion is met.

The active set is changed only when a free variable hits its bounds during an iteration or the optimality condition is met for the free variables but not for all variables in IA, the active set. In the latter case, a variable that violates the optimality condition will be dropped out of IA. For more details on the modified Newton method and line search, see Dennis and Schnabel (1983). For more detailed information on active set strategy, see Gill and Murray (1976).

Since a finite-difference method is used to estimate the Hessian for some single precision calculations, an inaccurate estimate of the Hessian may cause the algorithm to terminate at a noncritical point. In such cases, high precision arithmetic is recommended. Also, whenever the exact Hessian can be easily provided, routine BCOAH should be used instead.

Comments

1. Workspace may be explicitly provided, if desired, by use of B2ODH/DB2ODH. The reference is:

CALL B2ODH (FCN, GRAD, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FVALUE, WK, IWK)

The additional arguments are as follows:

WK — Real work vector of length N * (N + 8). WK contains the following information on output: The second N locations contain the last step taken. The third N locations contain the last Newton step. The fourth N locations contain an estimate of the gradient at the solution. The final N2 locations contain the Hessian at the approximate solution.

IWK — Integer work vector of length N.

2. Informational errors

Type | Code | Description |
---|---|---|

3 | 1 | Both the actual and predicted relative reductions in the function are less than or equal to the relative function convergence tolerance. |

4 | 2 | The iterates appear to be converging to a noncritical point. |

4 | 3 | Maximum number of iterations exceeded. |

4 | 4 | Maximum number of function evaluations exceeded. |

4 | 5 | Maximum number of gradient evaluations exceeded. |

4 | 6 | Five consecutive steps have been taken with the maximum step length. |

2 | 7 | Scaled step tolerance satisfied; the current point may be an approximate local solution, or the algorithm is making very slow progress and is not near a solution, or STEPTL is too big. |

4 | 7 | Maximum number of Hessian evaluations exceeded. |

3. The first stopping criterion for BCODH occurs when the norm of the gradient is less than the given gradient tolerance (RPARAM(1)). The second stopping criterion for BCODH occurs when the scaled distance between the last two steps is less than the step tolerance (RPARAM(2)).

4. If the default parameters are desired for BCODH, then set IPARAM(1) to zero and call the routine BCODH. Otherwise, if any nondefault parameters are desired for IPARAM or RPARAM; then the following steps should be taken before calling BCODH:

CALL U4INF (IPARAM, RPARAM)

Set nondefault values for desired IPARAM, RPARAM elements.

Note that the call to U4INF will set IPARAM and RPARAM to their default values so only nondefault values need to be set above. |

The following is a list of the parameters and the default values:

IPARAM — Integer vector of length 7.

IPARAM(1) = Initialization flag.

IPARAM(2) = Number of good digits in the function.

Default: Machine dependent.

Default: Machine dependent.

IPARAM(3) = Maximum number of iterations.

Default: 100.

Default: 100.

IPARAM(4) = Maximum number of function evaluations.

Default: 400.

Default: 400.

IPARAM(5) = Maximum number of gradient evaluations.

Default: 400.

Default: 400.

IPARAM(6) = Hessian initialization parameter.

Default: Not used in BCODH.

Default: Not used in BCODH.

IPARAM(7) = Maximum number of Hessian evaluations.

Default: 100.

Default: 100.

RPARAM — Real vector of length 7.

RPARAM(1) = Scaled gradient tolerance.

The i-th component of the scaled gradient at x is calculated as

The i-th component of the scaled gradient at x is calculated as

where g = ∇f (x), s = XSCALE, and fs = FSCALE.

Default:

Default:

in double where ɛ is the machine precision.

RPARAM(2) = Scaled step tolerance. (STEPTL)

The i-th component of the scaled step between two points x and y is computed as

The i-th component of the scaled step between two points x and y is computed as

where s = XSCALE.

Default: ɛ2/3 where ɛ is the machine precision.

Default: ɛ2/3 where ɛ is the machine precision.

RPARAM(3) = Relative function tolerance.

Default: max(10-10, ɛ2/3), max (10-20, ɛ2/3) in double where ɛ is the machine precision.

Default: max(10-10, ɛ2/3), max (10-20, ɛ2/3) in double where ɛ is the machine precision.

RPARAM(4) = Absolute function tolerance.

Default: Not used in BCODH.

Default: Not used in BCODH.

RPARAM(5) = False convergence tolerance.

Default: 100ɛ where ɛ is the machine precision.

Default: 100ɛ where ɛ is the machine precision.

RPARAM(6) = Maximum allowable step size.

Default: 1000 max(ɛ1, ɛ2) where

Default: 1000 max(ɛ1, ɛ2) where

ɛ2 = ∥s∥2, s = XSCALE, and t = XGUESS.

RPARAM(7) = Size of initial trust region radius.

Default: based on the initial scaled Cauchy step.

Default: based on the initial scaled Cauchy step.

If double precision is required, then DU4INF is called and RPARAM is declared double precision.

5. Users wishing to override the default print/stop attributes associated with error messages issued by this routine are referred to “Error Handling” in the Introduction.

Example

The problem

is solved with an initial guess (-1.2, 1.0), and default values for parameters.

USE BCODH_INT

USE UMACH_INT

IMPLICIT NONE

INTEGER N

PARAMETER (N=2)

!

INTEGER IP, IPARAM(7), L, NOUT

REAL F, X(N), XGUESS(N), XLB(N), XUB(N)

EXTERNAL ROSBRK, ROSGRD

!

DATA XGUESS/-1.2E0, 1.0E0/

DATA XLB/-2.0E0, -1.0E0/, XUB/0.5E0, 2.0E0/

!

IPARAM(1) = 0

IP = 0

! Minimize Rosenbrock function using

! initial guesses of -1.2 and 1.0

CALL BCODH (ROSBRK, ROSGRD, IP, XLB, XUB, X, XGUESS=XGUESS, &

IPARAM=IPARAM, FVALUE=F)

! Print results

CALL UMACH (2, NOUT)

WRITE (NOUT,99999) X, F, (IPARAM(L),L=3,5)

!

99999 FORMAT (' The solution is ', 6X, 2F8.3, //, ' The function ', &

'value is ', F8.3, //, ' The number of iterations is ', &

10X, I3, /, ' The number of function evaluations is ', &

I3, /, ' The number of gradient evaluations is ', I3)

!

END

!

SUBROUTINE ROSBRK (N, X, F)

INTEGER N

REAL X(N), F

!

F = 1.0E2*(X(2)-X(1)*X(1))**2 + (1.0E0-X(1))**2

!

RETURN

END

SUBROUTINE ROSGRD (N, X, G)

INTEGER N

REAL X(N), G(N)

!

G(1) = -4.0E2*(X(2)-X(1)*X(1))*X(1) - 2.0E0*(1.0E0-X(1))

G(2) = 2.0E2*(X(2)-X(1)*X(1))

!

RETURN

END

Output

The solution is 0.500 0.250

The function value is 0.250

The number of iterations is 17

The number of function evaluations is 26

The number of gradient evaluations is 18