An algorithm for general-purpose unconstrained non-linear optimization. The algorithm is of quasi-Newton type with BFGS updating of the inverse Hessian and soft line search with a trust region type monitoring of the input to the line search algorithm. The interface of ‘ucminf’ is designed for easy interchange with ‘optim’.
Usage
ucminf(par, fn, gr = NULL, ..., control = list(), hessian = 0)Arguments
- par
Initial estimate of minimum for
fn.- fn
Objective function to be minimized.
- gr
Gradient of objective function. If
NULLa finite difference approximation is used.- ...
Optional arguments passed to the objective and gradient functions.
- control
A list of control parameters. See ‘Details’.
- hessian
Integer value:
- 1
Returns a numerical approximation of the Hessian using ‘hessian’ in the package ‘numDeriv’.
- 2
Returns final approximation of the inverse Hessian based on the series of BFGS updates during optimization.
- 3
Same at 2, but will also return the Hessian (the inverse of 2).
If a
TRUEorFALSEvalue is given it will switch between option 1 or 0.
Value
ucminf returns a list of class 'ucminf'
containing the following elements:
- par
Computed minimizer.
- value
Objective function value at computed minimizer.
- convergence
Flag for reason of termination:
- 1
Stopped by small gradient (grtol).
- 2
Stopped by small step (xtol).
- 3
Stopped by function evaluation limit (maxeval).
- 4
Stopped by zero step from line search
- -2
Computation did not start: length(par) = 0.
- -4
Computation did not start: stepmax is too small.
- -5
Computation did not start: grtol or xtol <= 0.
- -6
Computation did not start: maxeval <= 0.
- -7
Computation did not start: given Hessian not pos. definite.
- message
String with reason of termination.
- hessian, invhessian
Estimate of (inv.) Hessian at computed minimizer. The type of estimate is given by the input argument ‘hessian’.
- invhessian.lt
The lower triangle of the final approximation to the inverse Hessian based on the series of BFGS updates during optimization.
- info
Information about the search:
- maxgradient
\(||F'(x)||_\infty\), the largest element in the absolute value of the gradient at the computed minimizer.
- laststep
Length of last step.
- stepmax
Final maximal allowed step length.
- neval
Number of calls to both objective and gradient function.
Details
The algorithm is documented in (Nielsen, 2000) (see References below) together with a comparison to the Fortran subroutine ‘MINF’ and the Matlab function ‘fminunc’. The implementation of ‘ucminf’ in R uses the original Fortran version of the algorithm.
The interface in R is designed so that it is very easy to switch
between using ‘ucminf’ and ‘optim’. The
arguments par, fn, gr, and hessian
are all the same (with a few extra options for hessian in
‘ucminf’). The difference is that there is no method
argument in ‘ucminf’ and that some of the components in the
control argument are different due to differences in the algorithms.
The algorithm can be given an initial estimate of the Hessian for the optimization and it is possible to get the final approximation of the Hessian based on the series of BFGS updates. This extra functionality may be useful for optimization in a series of related problems.
The functions fn and gr can return Inf or NaN
if the functions cannot be evaluated at the supplied value, but the
functions must be computable at the initial value. The functions
are not allowed to return NA. Any names given to par will be
copied to the vectors passed to fn and gr. No
other attributes of par are copied over.
The control argument is a list that can supply any of the
following components:
grtolThe algorithm stops when \(||F'(x)||_\infty \leq \) grtol, that is when the largest absolute value of the gradient is less than grtol. Default value is
grtol = 1e-6.xtolThe algorithm stops when \(||x-x_p||_2 \leq \textrm{xtol}\cdot(\textrm{xtol} + ||x||_2)\), where \(x_p\) and \(x\) are the previous and current estimate of the minimizer. Thus the algorithm stops when the last relative step length is sufficiently small. Default value is
xtol = 1e-12.stepmaxInitial maximal allowed step length (radius of trust-region). The value is updated during the optimization. Default value is
stepmax = 1.maxevalThe maximum number of function evaluations. A function evaluation is counted as one evaluation of the objective function and of the gradient function. Default value is
maxeval = 500.gradEither ‘forward’ or ‘central’. Controls the type of finite difference approximation to be used for the gradient if no gradient function is given in the input argument ‘gr’. Default value is
grad = 'forward'.gradstepVector of length 2. The step length in finite difference approximation for the gradient. Step length is \(|x_i|\cdot\textrm{gradstep[1]+gradstep[2]}\). Default value is
gradstep = c(1e-6, 1e-8).invhessian.ltA vector with an initial approximation to the lower triangle of the inverse Hessian. If not given, the inverse Hessian is initialized as the identity matrix. If
H0is the initial hessian matrix then the lower triangle of the inverse ofH0can be found asinvhessian.lt = solve(H0)[lower.tri(H0,diag=TRUE)].
References
Nielsen, H. B. (2000) ‘UCMINF - An Algorithm For Unconstrained, Nonlinear Optimization’, Report IMM-REP-2000-19, Department of Mathematical Modelling, Technical University of Denmark. http://www.imm.dtu.dk/documents/ftp/tr00/tr19_00.pdf
The original Fortran source code was found at
http://www2.imm.dtu.dk/projects/hbn_software/ucminf.f.
(That URL is no longer available but archived at
https://web.archive.org/web/20050418082240/http://www.imm.dtu.dk/~hbn/Software/ucminf.f
– Dr Nielsen passed away in 2015). The code has been slightly modified in
this package to be suitable for use with R.
The general structure of the implementation in R is based on the package ‘FortranCallsR’ by Diethelm Wuertz.
Author
‘UCMINF’ algorithm design and Fortran code by Hans Bruun Nielsen.
K Hervé Dakpo took over maintenance of the package in May. 2023.
Implementation in R by Stig B. Mortensen, stigbm@gmail.com.
Modifications by Douglas Bates bates@stat.wisc.edu, Nov. 2010, to support nested optimization and correct issues with printing on Windows.
Examples
## Rosenbrock Banana function
fR <- function(x) (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
gR <- function(x) c(-400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]),
200 * (x[2] - x[1] * x[1]))
# Find minimum and show trace
ucminf(par = c(2,.5), fn = fR, gr = gR, control = list(trace = 1))
#> neval = 1, F(x) = 1.2260e+03, max|g(x)| = 2.8020e+03
#> x = 2.0000e+00, 5.0000e-01
#> Line search: alpha = 1.0000e+00, dphi(0) =-2.8881e+03, dphi(1) =-1.4263e+02
#> neval = 2, F(x) = 1.0123e+01, max|g(x)| = 1.3111e+02
#> x = 1.0298e+00, 7.4237e-01
#> Line search: alpha = 1.0000e+00, dphi(0) =-3.1743e+01, dphi(1) = 1.0180e+01
#> neval = 3, F(x) = 1.7049e+00, max|g(x)| = 6.3969e+01
#> x = 1.2600e+00, 1.7155e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-2.5788e+00, dphi(1) =-5.6182e-01
#> neval = 4, F(x) = 1.1612e-01, max|g(x)| = 1.2343e+01
#> x = 1.2174e+00, 1.5083e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-1.5867e-01, dphi(1) = 1.2108e-02
#> neval = 5, F(x) = 4.2253e-02, max|g(x)| = 1.8638e+00
#> x = 1.2033e+00, 1.4449e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-1.1826e-03, dphi(1) =-3.2371e-04
#> neval = 6, F(x) = 4.1500e-02, max|g(x)| = 8.6681e-01
#> x = 1.2035e+00, 1.4474e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-5.9673e-04, dphi(1) =-4.7194e-04
#> neval = 7, F(x) = 4.0965e-02, max|g(x)| = 4.8839e-01
#> x = 1.2024e+00, 1.4456e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-3.9731e-03, dphi(1) =-2.3018e-03
#> neval = 8, F(x) = 3.7853e-02, max|g(x)| = 8.5215e-01
#> x = 1.1928e+00, 1.4254e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-8.0453e-03, dphi(1) =-6.3954e-03
#> neval = 9, F(x) = 3.0800e-02, max|g(x)| = 2.0990e+00
#> x = 1.1676e+00, 1.3685e+00
#> Line search: alpha = 8.2084e-01, dphi(0) =-4.4175e-02, dphi(1) = 1.8746e-02
#> neval = 11, F(x) = 4.8486e-03, max|g(x)| = 2.2862e+00
#> x = 1.0458e+00, 1.0884e+00
#> Line search: alpha = 3.8293e-01, dphi(0) =-4.8734e-03, dphi(1) = 4.6817e-04
#> neval = 13, F(x) = 4.0485e-03, max|g(x)| = 1.1863e+00
#> x = 1.0584e+00, 1.1177e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-6.4354e-04, dphi(1) =-5.6879e-04
#> neval = 14, F(x) = 3.4426e-03, max|g(x)| = 1.1238e+00
#> x = 1.0535e+00, 1.1074e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-4.7371e-03, dphi(1) =-1.0920e-03
#> neval = 15, F(x) = 6.1678e-04, max|g(x)| = 7.3075e-01
#> x = 1.0180e+00, 1.0347e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-7.9043e-04, dphi(1) =-2.5377e-04
#> neval = 16, F(x) = 1.0437e-04, max|g(x)| = 1.6394e-01
#> x = 1.0096e+00, 1.0189e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-1.8089e-04, dphi(1) =-1.8237e-05
#> neval = 17, F(x) = 5.8219e-06, max|g(x)| = 9.1455e-02
#> x = 1.0009e+00, 1.0016e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-1.3102e-05, dphi(1) = 2.0222e-06
#> neval = 18, F(x) = 2.9162e-07, max|g(x)| = 1.7185e-02
#> x = 1.0003e+00, 1.0007e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-5.9332e-07, dphi(1) = 1.1234e-08
#> neval = 19, F(x) = 1.2578e-10, max|g(x)| = 2.0751e-04
#> x = 9.9999e-01, 9.9998e-01
#> Line search: alpha = 1.0000e+00, dphi(0) =-2.5270e-10, dphi(1) = 1.1297e-12
#> neval = 20, F(x) = 3.5670e-15, max|g(x)| = 2.0836e-06
#> x = 1.0000e+00, 1.0000e+00
#> Line search: alpha = 1.0000e+00, dphi(0) =-7.1150e-15, dphi(1) =-1.8980e-17
#> Optimization has converged
#> Stopped by small gradient (grtol).
#> maxgradient laststep stepmax neval
#> 1.020598e-08 6.480989e-08 1.225000e-01 2.100000e+01
#> $par
#> [1] 1 1
#>
#> $value
#> [1] 8.073499e-20
#>
#> $convergence
#> [1] 1
#>
#> $message
#> [1] "Stopped by small gradient (grtol)."
#>
#> $invhessian.lt
#> [1] 0.4975315 0.9951877 1.9956183
#>
#> $info
#> maxgradient laststep stepmax neval
#> 1.020598e-08 6.480989e-08 1.225000e-01 2.100000e+01
#>
#> attr(,"class")
#> [1] "ucminf"