Skip to content

Commit

Permalink
multistart updates + control parameter tuning
Browse files Browse the repository at this point in the history
  • Loading branch information
@JorisChau
JorisChau committed Apr 16, 2024
1 parent 17ef629 commit 0442f0c
Show file tree
Hide file tree
Showing 5 changed files with 383 additions and 150 deletions.
62 changes: 40 additions & 22 deletions R/nls.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -383,15 +383,15 @@ gsl_nls.formula <- function(fn, data = parent.frame(), start,
## Handle infinite/missing start values
if(is.vector(start)) {
if(any((napars <- sapply(start, Negate(is.finite))))) {
start <- vapply(names(start), function(nm) if(napars[[nm]]) c(0, 1) else rep(start[[nm]], 2L), numeric(2))
start <- vapply(names(start), function(nm) if(napars[[nm]]) c(-0.1, 0.75) else rep(start[[nm]], 2L), numeric(2))
.has_start <- sapply(!napars, rep, times = 2L)
} else {
.has_start <- !napars ## not used (single-start)
}
} else if(is.matrix(start)) {
if(any(napars <- apply(start, 2L, Negate(is.finite)))) {
start[1L, napars[1L, ]] <- 0
start[2L, napars[2L, ]] <- 1
start[1L, napars[1L, ]] <- -0.1
start[2L, napars[2L, ]] <- 0.75
}
.has_start <- !napars
} else {
Expand Down Expand Up @@ -501,9 +501,18 @@ gsl_nls.formula <- function(fn, data = parent.frame(), start,
.lupars[2, match(names(upper), pnames)] <- unname(upper)
if(any(.lupars[1L, ] > .lupars[2L, ]))
stop("Parameter lower bounds cannot be larger than upper bounds")
if(is.matrix(start) && (any(start[1, ] < .lupars[1, ] | start[2, ] > .lupars[2, ])))
stop("Starting parameter ranges must be contained within 'lower' and 'upper' bounds")
else if(!is.matrix(start) && any(vapply(startnames, function(nm) start[[nm]] < .lupars[1, nm] || start[[nm]] > .lupars[2, nm], logical(1))))
if(is.matrix(start)) {
if(any(!.has_start[1, ])) {
start1 <- start[1, !.has_start[1, ]]
start[1, !.has_start[1, ]] <- pmax(start1, .lupars[1, !.has_start[1, ]])
start[2, !.has_start[2, ]] <- start[2, !.has_start[2, ]] + (start[1, !.has_start[1, ]] - start1)
}
if(any(!.has_start[2, ])) {
start[2, !.has_start[2, ]] <- pmin(start[2, !.has_start[2, ]], .lupars[2, !.has_start[2, ]])
}
if(any(start[1, ] < .lupars[1, ] | start[2, ] > .lupars[2, ]))
stop("Starting parameter ranges must be contained within 'lower' and 'upper' bounds")
} else if(!is.matrix(start) && any(vapply(startnames, function(nm) start[[nm]] < .lupars[1, nm] || start[[nm]] > .lupars[2, nm], logical(1))))
stop("Starting parameters must be contained within 'lower' and/or 'upper' bounds")
if(all(is.infinite(.lupars))) {
.lupars <- NULL
Expand Down Expand Up @@ -616,7 +625,7 @@ gsl_nls.formula <- function(fn, data = parent.frame(), start,
.ctrl$solver <- match.arg(.ctrl$solver, c("qr", "cholesky", "svd"))
.ctrl$fdtype <- match.arg(.ctrl$fdtype, c("forward", "center"))
stopifnot(
is.numeric(.ctrl$maxiter), length(.ctrl$maxiter) == 1, .ctrl$maxiter > 0,
is.numeric(.ctrl$maxiter), length(.ctrl$maxiter) == 1, .ctrl$maxiter >= 1,
is.numeric(.ctrl$factor_up), length(.ctrl$factor_up) == 1, .ctrl$factor_up > 0,
is.numeric(.ctrl$factor_down), length(.ctrl$factor_down) == 1, .ctrl$factor_down > 0,
is.numeric(.ctrl$avmax), length(.ctrl$avmax) == 1, .ctrl$avmax > 0,
Expand Down Expand Up @@ -734,15 +743,15 @@ gsl_nls.function <- function(fn, y, start,
## Handle infinite/missing start values
if(is.vector(.start)) {
if(any((napars <- sapply(.start, Negate(is.finite))))) {
.start <- vapply(names(.start), function(nm) if(napars[[nm]]) c(0, 1) else rep(.start[[nm]], 2L), numeric(2))
.start <- vapply(names(.start), function(nm) if(napars[[nm]]) c(-0.1, 0.75) else rep(.start[[nm]], 2L), numeric(2))
.has_start <- sapply(!napars, rep, times = 2L)
} else {
.has_start <- !napars ## not used (single-start)
}
} else if(is.matrix(.start)) {
if(any(napars <- apply(.start, 2L, Negate(is.finite)))) {
.start[1L, napars[1L, ]] <- 0
.start[2L, napars[2L, ]] <- 1
.start[1L, napars[1L, ]] <- -0.1
.start[2L, napars[2L, ]] <- 0.75
}
.has_start <- !napars
} else {
Expand Down Expand Up @@ -799,9 +808,18 @@ gsl_nls.function <- function(fn, y, start,
.lupars[2, match(names(upper), pnames)] <- unname(upper)
if(any(.lupars[1L, ] > .lupars[2L, ]))
stop("Parameter lower bounds cannot be larger than upper bounds")
if(is.matrix(.start) && (any(.start[1, ] < .lupars[1, ] | .start[2, ] > .lupars[2, ])))
stop("Starting parameter ranges must be contained within 'lower' and 'upper' bounds")
else if(!is.matrix(.start) && any(vapply(pnames, function(nm) .start[[nm]] < .lupars[1, nm] || .start[[nm]] > .lupars[2, nm], logical(1))))
if(is.matrix(.start)) {
if(any(!.has_start[1, ])) {
start1 <- .start[1, !.has_start[1, ]]
.start[1, !.has_start[1, ]] <- pmax(start1, .lupars[1, !.has_start[1, ]])
.start[2, !.has_start[2, ]] <- .start[2, !.has_start[2, ]] + (.start[1, !.has_start[1, ]] - start1)
}
if(any(!.has_start[2, ])) {
.start[2, !.has_start[2, ]] <- pmin(.start[2, !.has_start[2, ]], .lupars[2, !.has_start[2, ]])
}
if(any(.start[1, ] < .lupars[1, ] | .start[2, ] > .lupars[2, ]))
stop("Starting parameter ranges must be contained within 'lower' and 'upper' bounds")
} else if(!is.matrix(.start) && any(vapply(pnames, function(nm) .start[[nm]] < .lupars[1, nm] || .start[[nm]] > .lupars[2, nm], logical(1))))
stop("Starting parameters must be contained within 'lower' and/or 'upper' bounds")
if(all(is.infinite(.lupars))) {
.lupars <- NULL
Expand Down Expand Up @@ -888,7 +906,7 @@ gsl_nls.function <- function(fn, y, start,
.ctrl$solver <- match.arg(.ctrl$solver, c("qr", "cholesky", "svd"))
.ctrl$fdtype <- match.arg(.ctrl$fdtype, c("forward", "center"))
stopifnot(
is.numeric(.ctrl$maxiter), length(.ctrl$maxiter) == 1, .ctrl$maxiter > 0,
is.numeric(.ctrl$maxiter), length(.ctrl$maxiter) == 1, .ctrl$maxiter >= 1,
is.numeric(.ctrl$factor_up), length(.ctrl$factor_up) == 1, .ctrl$factor_up > 0,
is.numeric(.ctrl$factor_down), length(.ctrl$factor_down) == 1, .ctrl$factor_down > 0,
is.numeric(.ctrl$avmax), length(.ctrl$avmax) == 1, .ctrl$avmax > 0,
Expand Down Expand Up @@ -1025,10 +1043,10 @@ gsl_nls.function <- function(fn, y, start,
#' defaults to \code{sqrt(.Machine$double.eps)}.
#' @param gtol numeric value, termination occurs when the relative size of the gradient of the sum of squared residuals is \code{<= gtol},
#' indicating a local minimum, defaults to \code{.Machine$double.eps^(1/3)}
#' @param mstart_n positive integer, number of quasirandom points drawn in each major iteration, parameter \code{N} in Hickernell and Yuan (1997). Default is 25.
#' @param mstart_n positive integer, number of quasirandom points drawn in each major iteration, parameter \code{N} in Hickernell and Yuan (1997). Default is 30.
#' @param mstart_p positive integer, number of iterations of inexpensive local search to concentrate the sample, parameter \code{p} in Hickernell and Yuan (1997). Default is 5.
#' @param mstart_q positive integer, number of points retained in the concentrated sample, parameter \code{q} in Hickernell and Yuan (1997). Default is 5.
#' @param mstart_r positive integer, scaling factor of number of stationary points determining when the multi-start algorithm terminates, parameter \code{r} in Hickernell and Yuan (1997). Default is 3.
#' @param mstart_q positive integer, number of points retained in the concentrated sample, parameter \code{q} in Hickernell and Yuan (1997). Default is \code{mstart_n \%/\% 10}..
#' @param mstart_r positive integer, scaling factor of number of stationary points determining when the multi-start algorithm terminates, parameter \code{r} in Hickernell and Yuan (1997). Default is 4.
#' If the starting ranges for one or more parameters are unbounded and updated dynamically, \code{mstart_r} is multiplied by a factor 10 to avoid early termination.
#' @param mstart_s positive integer, minimum number of iterations a point needs to be retained before starting an efficient local search, parameter \code{s} in Hickernell and Yuan (1997). Default is 2.
#' @param mstart_tol numeric value, multiplicative tolerance \code{(1 + mstart_tol)} used as criterion to start an efficient local search (epsilon in Algorithm 2.1, Hickernell and Yuan (1997)).
Expand Down Expand Up @@ -1073,17 +1091,17 @@ gsl_nls.function <- function(fn, y, start,
gsl_nls_control <- function(maxiter = 100, scale = "more", solver = "qr",
fdtype = "forward", factor_up = 2, factor_down = 3, avmax = 0.75,
h_df = sqrt(.Machine$double.eps), h_fvv = 0.02, xtol = sqrt(.Machine$double.eps),
ftol = sqrt(.Machine$double.eps), gtol = .Machine$double.eps^(1/3),
mstart_n = 25, mstart_p = 5, mstart_q = 5, mstart_r = 3, mstart_s = 2,
mstart_tol = 0.5, mstart_maxiter = maxiter %/% 10,
mstart_maxstart = 1000, mstart_minsp = 1) {
ftol = sqrt(.Machine$double.eps), gtol = sqrt(.Machine$double.eps),
mstart_n = 30, mstart_p = 5, mstart_q = mstart_n %/% 10, mstart_r = 4, mstart_s = 2,
mstart_tol = 0.25, mstart_maxiter = maxiter %/% 10,
mstart_maxstart = 250, mstart_minsp = 1) {

scale <- match.arg(scale, c("more", "levenberg", "marquardt"))
solver <- match.arg(solver, c("qr", "cholesky", "svd"))
fdtype <- match.arg(fdtype, c("forward", "center"))

stopifnot(
is.numeric(maxiter), length(maxiter) == 1, maxiter > 0,
is.numeric(maxiter), length(maxiter) == 1, maxiter >= 1,
is.numeric(factor_up), length(factor_up) == 1, factor_up > 0,
is.numeric(factor_down), length(factor_down) == 1, factor_down > 0,
is.numeric(avmax), length(avmax) == 1, avmax > 0,
Expand Down
4 changes: 2 additions & 2 deletions R/nls_large.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -364,7 +364,7 @@ gsl_nls_large.formula <- function(fn, data = parent.frame(), start,
.ctrl$solver <- "cholesky" ## fixed
.ctrl$fdtype <- match.arg(.ctrl$fdtype, c("forward", "center"))
stopifnot(
is.numeric(.ctrl$maxiter), length(.ctrl$maxiter) == 1, .ctrl$maxiter > 0,
is.numeric(.ctrl$maxiter), length(.ctrl$maxiter) == 1, .ctrl$maxiter >= 1,
is.numeric(.ctrl$factor_up), length(.ctrl$factor_up) == 1, .ctrl$factor_up > 0,
is.numeric(.ctrl$factor_down), length(.ctrl$factor_down) == 1, .ctrl$factor_down > 0,
is.numeric(.ctrl$avmax), length(.ctrl$avmax) == 1, .ctrl$avmax > 0,
Expand Down Expand Up @@ -541,7 +541,7 @@ gsl_nls_large.function <- function(fn, y, start,
.ctrl$solver <- "cholesky" ## fixed
.ctrl$fdtype <- match.arg(.ctrl$fdtype, c("forward", "center"))
stopifnot(
is.numeric(.ctrl$maxiter), length(.ctrl$maxiter) == 1, .ctrl$maxiter > 0,
is.numeric(.ctrl$maxiter), length(.ctrl$maxiter) == 1, .ctrl$maxiter >= 1,
is.numeric(.ctrl$factor_up), length(.ctrl$factor_up) == 1, .ctrl$factor_up > 0,
is.numeric(.ctrl$factor_down), length(.ctrl$factor_down) == 1, .ctrl$factor_down > 0,
is.numeric(.ctrl$avmax), length(.ctrl$avmax) == 1, .ctrl$avmax > 0,
Expand Down
35 changes: 30 additions & 5 deletions src/gsl_nls.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
Expand Down Expand Up @@ -35,8 +36,11 @@ typedef struct
gsl_matrix *mx; // multistart initial value matrix
gsl_vector *mp; // multistart parameter placeholder vector
gsl_vector *mpopt; // multistart optimal parameter vector
gsl_vector *mpopt1; // multistart back-up parameter vector
gsl_vector *diag; // diagonal of scaling matrix D_k
gsl_matrix *lu; // lower/upper parameter bounds
gsl_matrix *JTJ; // workspace hessian matrix Jt * J
gsl_vector *workn; // workspace n-length vector
} pdata;

typedef struct
Expand All @@ -56,6 +60,15 @@ typedef struct
int startisnum; // are start values numeric?
} fdata;

/* store info found local optima */
typedef struct
{
double *mpall; // all found local optima
double *mpradii; // trust region radii of all found local optima
R_len_t mpcount; // current number of local optima
R_len_t mpmax; // maximum number of local optima
} mpoptinfo;

typedef struct
{
R_len_t n; // number of quasirandom samples
Expand All @@ -67,22 +80,27 @@ typedef struct
R_len_t minsp; // minimum number of stationary points
double r; // stopping criterion nwsp >= r * nsp
double tol; // start local search if mssr < (1 + tol) * mssropt
double dtol; // discard point if det(J^T * J) < dtol
int *ntix; // persistent index counter across major iters
double *qmp; // quasirandom samples
double *startptr; // pointer to start value matrix
double *start; // start value workspace
double *maxlims; // store maximum parameter limits
int *mssr_order; // ordered ssr's
int mstop; // stop status
R_len_t mstarts; // current number of major iteration
R_len_t nsp; // current number of stationary points
R_len_t nwsp; // current number of worse stationary points
double mssropt; // current optimal ssr
double ssrconv; // achieved convergence tolerance
Rboolean all_start; // flag if all fixed lower/upper sampling limits
int *has_start; // flags fixed lower/upper sampling limits
double mssropt[2]; // current optimal ssr + back-up
double ssrconv[2]; // current convergence tolerance + back-up
int all_start; // flag if any dynamic lower/upper sampling limits
int *has_start; // flags fixed lower/upper sampling limits
int *luchange; // force change in lower/upper sampling limits
double rejectscl; // scaling factor parameter rejection limits
mpoptinfo *mpopt; // store all found local optima (currently not)
} mdata;



/* from multifit_nlinear/trust.c */
typedef struct
{
Expand Down Expand Up @@ -115,6 +133,13 @@ int trust_iterate_lu(void *vstate,
gsl_vector *dx,
const gsl_matrix *lu);

double det_cholesky_jtj(gsl_matrix *J, gsl_matrix *JTJ);

double det_eval_jtj(const gsl_multifit_nlinear_parameters params,
const gsl_vector *swts, gsl_multifit_nlinear_fdf *fdf,
const gsl_vector *x, gsl_vector *f, gsl_matrix *J,
gsl_matrix *JTJ, gsl_vector *workn);

SEXP C_nls(SEXP fn, SEXP y, SEXP jac, SEXP fvv, SEXP env, SEXP start, SEXP swts, SEXP lupars, SEXP control_int, SEXP control_dbl, SEXP has_start);

SEXP C_nls_internal(void *data);
Expand Down
Loading

0 comments on commit 0442f0c

Please sign in to comment.