Run nonlinear minimizer using trust region algorithm with conjugate gradient search directions and quasi-Hessian updates.
| x | A numeric vector of starting values for the optimizer. |
|---|---|
| fn | An R function that takes |
| gr | An R function that takes x as its first argument. Returns a
numeric vector that is the gradient of |
| hs | An R function that takes x as its first argument.
It returns a Hessian matrix object of class |
| method | Valid arguments are |
| control | A list containing control parameters for the optimizer. See details. |
| ... | Additional arguments passed to |
List containing the following items:
Value of the objective function
Parameter vector at the optimum
Gradient at the optimum
Estimate of the Hessian at the optimum (as class
symmetricMatrix, returned only for Sparse method).
Number of iterations before stopping
A message describing the last state of the iterator
For the Sparse method only, the number of nonzero elements in the lower triangle of the Hessian
.
The following sections explain how to use the package as a whole.
The control list should include the following parameters.
Initial radius of the trust region. Default is 5. If the algorithm returns non-finite values of the objective function early in the process, try a lower number.
Minimum radius of trust region. Algorithm will terminate if radius is below this value. This is because it may not be possible to get the norm of the gradient smaller than prec, and this is another way to get the algorithm to stop.
tolerance for the conjugate gradient algorithm that is used for the trust region subproblem. Set it to something very small. Default is sqrt(.Machine$double.eps)
Precision for how close the norm of the gradient at the solution should be to zero, before the algorithm halts. It is possible that the algorithm will not get that far, so it will also stop when the radius of the trust region is smaller thanstop.trust.radius. If the trust region radius collapses, but the norm of the gradient really isn't close to zero, then something terrible has happened.
An integer. The frequency at which the algorithm
will display the current iteration number or function value, among
other things (see report.level). Defaults to 1.
The amount of detail in each report. Defaults to 2.
The number of significant digits used in each report. Defaults to 5.
The number of lines of iterations before the report column headers are reprinted. Defaults to 25.
Maximum number of iterations. Defaults to 100.
When the algorithm decides to shrink the trust region, it will multiply the trust radius by this factor. Defaults to 0.5.
When the algorithm decides to expand the trust region, it will multiply the algorithm by this factor. Defaults to 3.
The algorithm with accept a proposed move if the ratio of the actual improvement in the objective function, to the predicted improvement from the trust region subproblem, is greater than this amount. Otherwise, the trust region will contract. Default is 0.25.
First criterion to determine if the trust region should expand. If the ratio of the actual and proposed improvements in the objective function is less than this factor, the algorithm will consider expanding the trust region. See expand.threshold.radius. Default is 0.8.
If the ratio of the actual and proposed improvement in the objective function is less than expand.threshold.ap, then, if the normed distance of the proposed move is greater than expand.threshold.radius, times the current trust region radius, the trust region will expand. Default is 0.8.
The algorithm will minimize fn times this
factor. If you want to maximize fn, this value should be negative
(usually -1). Default is 1.
Frequency at which the preconditioner for the conjugate gradiate estimation of the trust region subproblem is reestimated. Preconditioners can help the convergence properties of the algorithm. Default is 1.
ID for choice of preconditioner. 0 is the
identity matrix (default), For the Sparse method, 1 is a modified Cholesky preconditioner. For the BFGS method, 1 is the full Cholesky decomposition. If you select 1 for the SR1 method, the algorithm will use the identity preconditioner instead.
Maximum number of conjugate gradient iterations to run when solving the trust region subproblem. A higher number will lead to more accurate solutions to the subproblem, but may also lead to longer run times. Defaults to 2000.
The report.level control parameter determines how much information is displayed each time the algorithm reports the current state. Possible values are
No information (a quiet run)
Current iteration number, and current value of the objective function.
Information from level 1, plus the current norm of the gradient and a status message.
Information from levels 1 and 2, plus the current normed radius of the trust region.
Information from levels 1, 2, and 3, plus information from each estimate of the trust region subproblem (number of conjugate gradient iterations and how/why the CG algorithm terminated).
Default level is 2. Levels 3 and 4 are available primarily for debugging purposes.
The algorithm will stop when one of the following three conditions are met:
The norm of the gradient, divided by the square root of the number of parameters, is less than prec.
The trust region collapse to a radius smaller than machine precision
The algorithm proposes zero or negative improvement in the objective function (should never happen)
The number of iterations reaches the control parameter maxit
If the algorithm appears to have stopped prematurely (i.e., the norm of the gradient is still too large), then one might just restart the algorithm. For the quasi-Newton algorithms (SR1 and BFGS), this will refresh the Hessian, and might allow more progress to be made.
Sometimes estimating the Hessian is easy (e.g., you have an analytic representation, or you are using some kind of algorithmic differentiation software). If you do not know the Hessian, but you do know the sparsity structure, try the sparseHessianFD package. The routines in sparseHessianFD compute the Hessian using finite differencing, but in a way that exploits the sparsity structure. In many cases, this can be faster than constructing an analytic Hessian for a large problem (e.g., when the Hessian has a block-arrow structure with a large number of blocks).
To use the sparseHessianFD package, you need to provide the row and column indices of the non-zero elements of the lower triangle of the Hessian. This structure cannot change during the course of the trust.optim routine. Also, you really should provide an analytic gradient. sparseHessianFD computes finite differences of the gradient, so if the gradient itself is finite-differenced, so much error is propagated through that the Hessians are nearly worthless close to the optimum.
Of course, sparseHessianFD is useful only for the Sparse method. That said, one may still get decent performance using these routines even if the Hessian is sparse, if the problem is not too large. Just treat the Hessian as if it were sparse.
if (FALSE) {
data(binary)
N <- length(binary$Y)
k <- NROW(binary$X)
start <- rep(0,(N+1)*k)
priors <- list(inv.Sigma = diag(k), inv.Omega = diag(k))
opt <- trust.optim(start, fn=binary.f,
gr = binary.grad,
hs = binary.hess,
method = "Sparse",
control = list(
report.precision=1L,
function.scale.factor=-1
),
data=binary, priors=priors
)
}