Besides trust(),
RcppTrust ships three more trust-region optimizers. Each
one has the same two layers as trust(): an R function for
interactive use, and a thread-safe C core with no R API calls and no
global state. The C core can run as an inner or outer optimizer inside
parallel C++ code such as nlmixr2est.
| R function | Ported from | Derivatives needed | Bounds | Thread-safe C entry point |
|---|---|---|---|---|
trust() |
CRAN trust (Geyer) |
gradient + Hessian | no | trust_solve_c() |
steihaug() |
Rust crate basin | gradient + Hessian, or Hessian-vector products | no | steihaug_solve_c() |
bobyqa() |
Powell’s BOBYQA Fortran, via CRAN minqa |
none | box | bobyqa_solve_c() |
newuoa() |
Powell’s NEWUOA Fortran, via CRAN minqa |
none | no | newuoa_solve_c() |
All four repeat the same loop: build a quadratic model of the objective, minimize it within a radius around the current point, compare the actual decrease with the predicted one, and grow or shrink the radius. They differ in where the model comes from and how the subproblem is solved.
trust() uses the exact gradient and
Hessian, and solves the subproblem nearly exactly with an
eigendecomposition. That is robust, but each iteration costs
O(n^3).steihaug() uses the same
derivative-based model, but solves the subproblem approximately by
conjugate gradients (CG). CG stops early when it reaches the
trust-region boundary or finds a direction of negative curvature. The
Hessian is never factorized. If you supply only Hessian-vector products
(hessvec), the Hessian is never formed at all.bobyqa() and newuoa()
never use derivatives. They fit the quadratic model through
npt previously evaluated points (usually
2n + 1), and take up the model’s remaining freedom with the
smallest change to its Hessian. The step is found by truncated CG, so
NEWUOA’s subproblem solver is close in spirit to Steihaug’s; BOBYQA’s
also respects box bounds. They track two radii: rho, a
lower bound that only shrinks, from rhobeg to
rhoend, and controls termination; and delta,
the step radius, which can grow and shrink above it.steihaug() takes separate objective, gradient and
Hessian functions, because basin requests each one separately. It
reproduces basin’s evaluation pattern:
fn and gr once at the start;hess once per outer iteration;fn at each trial point;gr after each accepted step.library(RcppTrust)
fr <- function(x) 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
grr <- function(x) c(-400 * x[1] * (x[2] - x[1]^2) - 2 * (1 - x[1]),
200 * (x[2] - x[1]^2))
hr <- function(x) matrix(c(1200 * x[1]^2 - 400 * x[2] + 2, -400 * x[1],
-400 * x[1], 200), 2, 2)
r <- steihaug(c(-1.2, 1), fr, grr, hr)
r[c("par", "value", "iterations", "message", "counts")]
#> $par
#> [1] 1 1
#>
#> $value
#> [1] 8.628062e-25
#>
#> $iterations
#> [1] 27
#>
#> $message
#> [1] "GradientTolerance"
#>
#> $counts
#> fn gr hess hessvec
#> 31 28 27 0In matrix-free mode you supply hessvec(x, v) instead of
hess. The iterates are the same, but CG calls
hessvec once per product instead of forming the matrix:
hv <- function(x, v) drop(hr(x) %*% v)
m <- steihaug(c(-1.2, 1), fr, grr, hessvec = hv)
identical(m$par, r$par)
#> [1] TRUE
m$counts
#> fn gr hess hessvec
#> 31 28 0 76control$trace = TRUE returns one row per subproblem
attempt, with the radius, predicted and actual reduction,
rho, whether the step was accepted, whether it hit the
boundary, and how many CG iterations it took:
tr <- steihaug(c(-1.2, 1), fr, grr, hr, control = list(trace = TRUE))$trace
head(tr)
#> iter radius value gradnorm trialValue preddiff rho stepnorm
#> 1 0 1.00 24.200000 232.867688 4.567782 18.0216125 1.0893708 0.15477985
#> 2 1 1.00 4.567782 30.944982 4.128383 0.4307893 1.0199868 0.02784227
#> 3 2 1.00 4.128383 1.948900 5.660014 1.8110241 -0.8457264 1.00000000
#> 4 2 0.25 4.128383 1.948900 3.679090 0.4389264 1.0236167 0.25000000
#> 5 3 0.50 3.679090 2.542429 3.023221 0.9038790 0.7256159 0.50000000
#> 6 4 0.50 3.023221 18.477217 2.717184 0.2977665 1.0277778 0.03223066
#> accept hitBoundary cgIter
#> 1 TRUE FALSE 1
#> 2 TRUE FALSE 1
#> 3 FALSE TRUE 2
#> 4 TRUE TRUE 2
#> 5 TRUE TRUE 2
#> 6 TRUE FALSE 1The C++ port reproduces basin’s arithmetic in the same order, with floating-point contraction turned off. On 33 test setups, each run in both exact and matrix-free mode, it matches a basin reference program bit for bit: every evaluation point, iteration count, termination reason and evaluation count.
The differences from basin are deliberate:
maxit is 100 (basin’s
executor allows 1000) and the absolute gradient tolerance is on at
1e-8 (off in basin). Every other tolerance is off unless
you set it in control, as in basin.bobyqa() and newuoa() are drop-in
replacements for minqa::bobyqa() and
minqa::newuoa(). The arguments, control
settings, defaults, warnings, printed output and return values are the
same:
fr <- function(x) 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
bobyqa(c(1, 2), fr, lower = c(0, 0), upper = c(4, 4))
#> parameter estimates: 0.999999968901681, 0.999999928305543
#> objective: 9.98796325533312e-15
#> number of function evaluations: 341
newuoa(c(-1.2, 1), fr)
#> parameter estimates: 0.999999953422979, 0.999999906557554
#> objective: 2.17773673529253e-15
#> number of function evaluations: 257Powell’s Fortran 77 code (the unmodified originals are archived in PRIMA),
in the modified form distributed with minqa, was translated
line by line into C++ (Powell, 2006; Powell, 2009). The port keeps
Powell’s variable names, workspace layout, goto control
flow and order of operations. Compared with minqa, the
parameters, objective values, evaluation counts and error codes are
bitwise identical:
extRosen <- function(x) {
n <- length(x)
sum(100 * (x[-1] - x[-n]^2)^2 + (1 - x[-n])^2)
}
ours <- bobyqa(rep(0.5, 6), extRosen, lower = -2, upper = 2,
control = list(npt = 13))
theirs <- minqa::bobyqa(rep(0.5, 6), extRosen, lower = -2, upper = 2,
control = list(npt = 13))
identical(unclass(ours), unclass(theirs))
#> [1] TRUEOne intentional difference. In minqa
1.2.8, BOBYQA’s RESCUE routine calls the objective as
CALFUN(N,X,IPRINT), but X is not defined
inside RESCUE. As a result, minqa evaluates
the objective at whatever is in uninitialized memory, and the same run
can give different results from one call to the next. The port evaluates
the point Powell’s original code intended, W(1..N).
RESCUE is only entered on badly conditioned problems, so
ordinary runs never hit it. The runs that do enter RESCUE
were checked against a copy of minqa with only that line
fixed, and matched it exactly.
UOBYQA from minqa is not included. It needs
(n+1)(n+2)/2 points, which grows as O(n^2), and NEWUOA
covers the same unconstrained, derivative-free case with far fewer
points.
Each solver has a plain C header, installed with the package:
minqa_types.h: bobyqa_solve_c(),
newuoa_solve_c(), minqa_result_free(),
minqa_options_default();steihaug_types.h: steihaug_solve_c(),
steihaug_result_free(),
steihaug_options_default().The conventions are the same as for trust_solve_c():
userdata pointer.*_result_free().For the minqa solvers, printing (iprint) goes through an
optional callback instead of Rprintf(), and is silent by
default.
A consuming package reaches these cores through the same header-only
pointer table described in vignette("RcppTrust"). After
iniRcppTrustPtrs() runs, bobyqa_solve_c_ptr,
newuoa_solve_c_ptr, minqa_result_free_ptr,
steihaug_solve_c_ptr and
steihaug_result_free_ptr are available alongside
trust_solve_c_ptr. New slots are only ever appended to the
table. When the installed RcppTrust is older than the
header a consumer was built against, the missing pointers are left
NULL, so check them before use.
The example below runs NEWUOA and matrix-free Steihaug on many starting points inside an OpenMP loop, through the resolved pointers:
cpp_code <- '
// [[Rcpp::depends(RcppTrust)]]
#include <Rcpp.h>
#include <vector>
extern "C" {
#define iniRcppTrustPtrs _vignette2_iniRcppTrustPtrs
#include <RcppTrust.h>
iniRcppTrust
}
// Derivative-free objective for NEWUOA.
extern "C" int rosen_f(int n, const double *x, double *f, void *) {
double t = x[1] - x[0] * x[0];
*f = 100.0 * t * t + (1.0 - x[0]) * (1.0 - x[0]);
return 0;
}
// Value/gradient on request, for Steihaug. The Hessian is never requested
// in matrix-free mode.
extern "C" int rosen_vg(int n, const double *x, double *value,
double *gradient, double *hessian, void *) {
double t = x[1] - x[0] * x[0];
if (value) *value = 100.0 * t * t + (1.0 - x[0]) * (1.0 - x[0]);
if (gradient) {
gradient[0] = -400.0 * x[0] * t - 2.0 * (1.0 - x[0]);
gradient[1] = 200.0 * t;
}
return 0;
}
extern "C" int rosen_hv(int n, const double *x, const double *v, double *hv,
void *) {
double h11 = 1200.0 * x[0] * x[0] - 400.0 * x[1] + 2.0, h12 = -400.0 * x[0];
hv[0] = h11 * v[0] + h12 * v[1];
hv[1] = h12 * v[0] + 200.0 * v[1];
return 0;
}
// [[Rcpp::export]]
Rcpp::NumericMatrix fit_many(SEXP ptrTable, int nStarts) {
_vignette2_iniRcppTrustPtrs(ptrTable); // normally done once in .onLoad()
Rcpp::NumericMatrix out(nStarts, 4);
std::vector<double> res(4 * nStarts);
// Build with OpenMP (e.g. -fopenmp in the package Makevars) to run
// this loop in parallel; without it the pragma is ignored.
#pragma omp parallel for num_threads(2)
for (int i = 0; i < nStarts; i++) {
double start[2] = {-1.2 + 0.05 * i, 1.0};
minqa_options_t mo = minqa_options_default(2, start);
minqa_result_t mr;
newuoa_solve_c_ptr(2, start, rosen_f, nullptr, &mo, &mr);
steihaug_options_t so = steihaug_options_default();
steihaug_result_t sr;
steihaug_solve_c_ptr(2, start, rosen_vg, rosen_hv, nullptr, &so, &sr);
res[4 * i] = mr.par[0];
res[4 * i + 1] = mr.feval;
res[4 * i + 2] = sr.argument[0];
res[4 * i + 3] = sr.iterations;
minqa_result_free_ptr(&mr);
steihaug_result_free_ptr(&sr);
}
for (int i = 0; i < nStarts; i++)
for (int j = 0; j < 4; j++) out(i, j) = res[4 * i + j];
Rcpp::colnames(out) = Rcpp::CharacterVector::create(
"newuoa_x1", "newuoa_feval", "steihaug_x1", "steihaug_iter");
return out;
}
'
Rcpp::sourceCpp(code = cpp_code)
head(fit_many(RcppTrust:::.RcppTrustPtr(), 16L))
#> newuoa_x1 newuoa_feval steihaug_x1 steihaug_iter
#> [1,] 0.9999998 257 1 27
#> [2,] 1.0000000 259 1 43
#> [3,] 0.9999999 262 1 36
#> [4,] 1.0000000 285 1 30
#> [5,] 0.9999996 282 1 29
#> [6,] 0.9999999 298 1 32Every call inside the loop gets its own options, result and workspace, so the threads need no coordination. The package’s test suite runs each solver on 64 starting points both sequentially and across OpenMP threads, and checks that the two sets of results are bitwise identical.
Powell, M. J. D. (2006). The NEWUOA software for unconstrained optimization without derivatives. In G. Di Pillo and M. Roma (eds), Large-Scale Nonlinear Optimization, 255-297. Springer. https://doi.org/10.1007/0-387-30065-1_16
Powell, M. J. D. (2009). The BOBYQA algorithm for bound constrained optimization without derivatives. Report DAMTP 2009/NA06, Centre for Mathematical Sciences, University of Cambridge.
Powell, M. J. D. Original Fortran 77 source code for NEWUOA and BOBYQA, archived at https://github.com/libprima/prima/tree/main/fortran/original.
Steihaug, T. (1983). The conjugate gradient method and trust regions in large scale optimization. SIAM Journal on Numerical Analysis, 20(3), 626-637. https://doi.org/10.1137/0720042