| Title: | Graphical and Numerical Checks for Mode-Finding Routines |
|---|---|
| Description: | Tools for checking that the output of an optimization algorithm is indeed at a local mode of the objective function. This is accomplished graphically by calculating all one-dimensional "projection plots" of the objective function, i.e., varying each input variable one at a time with all other elements of the potential solution being fixed. The numerical values in these plots can be readily extracted for the purpose of automated and systematic unit-testing of optimization routines. |
| Authors: | Martin Lysy [aut, cre] |
| Maintainer: | Martin Lysy <[email protected]> |
| License: | GPL-3 |
| Version: | 1.0.1 |
| Built: | 2024-10-24 05:44:43 UTC |
| Source: | https://github.com/mlysy/optimcheck |
Tools for checking that the output of an optimization algorithm is indeed at a local mode of the objective function. This is accomplished graphically by calculating all one-dimensional "projection plots" of the objective function, i.e., varying each input variable one at a time with all other elements of the potential solution being fixed. The numerical values in these plots can be readily extracted for the purpose of automated and systematic unit-testing of optimization routines.
Maintainer: Martin Lysy [email protected]
Useful links:
# example: logistic regression ilogit <- binomial()$linkinv # generate data p <- sample(2:10,1) # number of parameters n <- sample(1000:2000,1) # number of observations X <- matrix(rnorm(n*p),n,p) # design matrix beta0 <- rnorm(p, sd = .1) # true parameter values y <- rbinom(n, size = 1, prob = ilogit(X %*% beta0)) # response # fit logistic regression bhat <- coef(glm(y ~ X - 1, family = binomial)) # check convergence # likelihood function loglik <- function(beta, y, X) { sum(dbinom(y, size = 1, prob = ilogit(X %*% beta), log = TRUE)) } # projection plots bnames <- parse(text = paste0("beta[", 1:p, "]")) system.time({ oproj <- optim_proj(xsol = bhat, fun = function(beta) loglik(beta, y, X), xnames = bnames, xlab = "Coefficient", ylab = "Loglikelihood") }) # numerical summary oproj # see ?summary.optproj for more information # elementwise differences between potential and optimal solution diff(oproj) # same as summary(oproj)$xdiff # refit general purpose optimizer starting from bhat # often faster than optim_proj, but less stable system.time({ orefit <- optim_refit(xsol = bhat, fun = function(beta) loglik(beta, y, X)) }) orefit# example: logistic regression ilogit <- binomial()$linkinv # generate data p <- sample(2:10,1) # number of parameters n <- sample(1000:2000,1) # number of observations X <- matrix(rnorm(n*p),n,p) # design matrix beta0 <- rnorm(p, sd = .1) # true parameter values y <- rbinom(n, size = 1, prob = ilogit(X %*% beta0)) # response # fit logistic regression bhat <- coef(glm(y ~ X - 1, family = binomial)) # check convergence # likelihood function loglik <- function(beta, y, X) { sum(dbinom(y, size = 1, prob = ilogit(X %*% beta), log = TRUE)) } # projection plots bnames <- parse(text = paste0("beta[", 1:p, "]")) system.time({ oproj <- optim_proj(xsol = bhat, fun = function(beta) loglik(beta, y, X), xnames = bnames, xlab = "Coefficient", ylab = "Loglikelihood") }) # numerical summary oproj # see ?summary.optproj for more information # elementwise differences between potential and optimal solution diff(oproj) # same as summary(oproj)$xdiff # refit general purpose optimizer starting from bhat # often faster than optim_proj, but less stable system.time({ orefit <- optim_refit(xsol = bhat, fun = function(beta) loglik(beta, y, X)) }) orefit
Elementwise difference between potential and optimal solutions.
## S3 method for class 'optcheck' diff(x, ...) ## S3 method for class 'summary.optcheck' diff(x, ...)## S3 method for class 'optcheck' diff(x, ...) ## S3 method for class 'summary.optcheck' diff(x, ...)
x |
Object of class |
... |
Further arguments to be passed to or from other methods. |
This function is simply a wrapper to summary(x)$xdiff and x$xdiff, for optcheck and summary.optcheck objects respectively.
A two-column matrix consisting of the absolute and relative differences between the potential and optimal solutions (xsol and xopt).
Given the objective function of an optimization problem and a potential solution, calculates "projection plots" along each coordinate of the solution vector, with all other coordinates being fixed at the input values.
optim_proj( xsol, fun, maximize = TRUE, xrng = 0.1, npts = 100, plot = TRUE, ... )optim_proj( xsol, fun, maximize = TRUE, xrng = 0.1, npts = 100, plot = TRUE, ... )
xsol |
Potential solution vector of length |
fun |
Objective function to be maximized (or minimized), with first argument the length- |
maximize |
Logical, whether a maximum or a minimum of the objective function is sought. |
xrng |
Optional specification of the range of each projection plot. Can be: (i) a |
npts |
Number of points in each projection plot. |
plot |
Logical, whether or not to display the projection plots or just return their contents. |
... |
Further arguments to pass to the |
An object of class optproj inheriting from optcheck (returned invisibly if plot = TRUE, with elements:
xsolThe potential solution.
ysolThe value of fun(xsol).
maximizeLogical indicating whether the potential solution should maximize or minimize the objective function.
xprojAn npts x nx matrix where each column is the x-axis of the projection plot along the given component of theta.
yprojAn npts x nx matrix where each column is the y-axis of the corresponding projection plot.
plot, summary, print, and diff methods for projection plots are available; see plot.optproj(), summary.optproj(), print.optproj(), and diff.optproj().
If the potential solution is indeed a local optimum of the objective function, and if it is used to initialize a second optimization, then original and "refined" solutions ought to be close.
optim_refit(xsol, fun, maximize = TRUE, maxit = 5000, reltol = 1e-08, xopt)optim_refit(xsol, fun, maximize = TRUE, maxit = 5000, reltol = 1e-08, xopt)
xsol |
Potential solution vector of length |
fun |
Objective function to be maximized (or minimized), with first argument the length- |
maximize |
Logical, whether a maximum or a minimum of the objective function is sought. |
maxit |
Maximum number of iterations for |
reltol |
Relative tolerance for convergence of |
xopt |
Optional refit solution calculated externally from an optimization algorithm of choice (see Details). |
By default, a so-called refined op(t)imization (or refit) test is performed by running the default Nelder-Mead simplex method provided by stats::optim(), initialized by the potential solution xsol. Only a simplified interface to stats::optim()'s control parameters are provided here.
Alternatively, the refit test can be performed with any optimization algorithm of choice. This is done externally, with the refined solution passed to optim_refit() via the argument xopt.
An object of class optrefit inheriting from optcheck, with elements:
xsolThe potential solution.
ysolThe value of fun(xsol).
maximizeLogical indicating whether the potential solution should maximize or minimize the objective function.
xoptThe solution found by the general-purpose optimizer.
yoptThe function value at the optimal solution, i.e., fun(xopt).
summary, print, and diff for optrefit objects are available; see summary.optrefit(), print.optrefit(), and diff.optrefit().
Projection plots for optimization routines.
## S3 method for class 'optproj' plot(x, xnames, xind, equalize = FALSE, layout, xlab, ylab, ...)## S3 method for class 'optproj' plot(x, xnames, xind, equalize = FALSE, layout, xlab, ylab, ...)
x |
An |
xnames |
Optional vector of names for the plot titles. |
xind |
Integer or logical vector of indices indicating which projections should be plotted. Defaults to all projection plots. |
equalize |
If |
layout |
Optional vector giving the number of rows and columns in the plot layout. For |
xlab, ylab
|
Outer x-axis and y-axis labels. |
... |
Further arguments to be passed to or from other methods. |
A grid of projection plots, with vertical lines at the potential solution.
optcheck and summary.optcheck objects.Print method for optcheck and summary.optcheck objects.
## S3 method for class 'optcheck' print(x, digits = max(3L, getOption("digits") - 3L), n = 5L, ...) ## S3 method for class 'summary.optcheck' print(x, digits = max(3L, getOption("digits") - 3L), n = 5L, ...)## S3 method for class 'optcheck' print(x, digits = max(3L, getOption("digits") - 3L), n = 5L, ...) ## S3 method for class 'summary.optcheck' print(x, digits = max(3L, getOption("digits") - 3L), n = 5L, ...)
x |
Object of class |
digits |
Number of digits to display. |
n |
Number of elements of solution vector to display (see Details). |
... |
Further arguments to be passed to or from other methods. |
The print methods for optcheck and summary.optcheck objects both display three-column matrix, consisting of the potential solution (xsol), the absolute difference between it and the optimal solution (xopt) return by either optim_proj() and optim_refit(), and the relative difference (R = (xopt - xsol)/|xsol|). Only the elemnts corresponding to the top-n relative differences are displayed.
Invisibly x itself.
summary method for projection plots.summary method for projection plots.
## S3 method for class 'optproj' summary(object, xnames, ...)## S3 method for class 'optproj' summary(object, xnames, ...)
object |
An |
xnames |
Optional vector of names for the elements of the potential solution. |
... |
Further arguments to be passed to or from other methods. |
The print methods for summary.optproj and optproj objects themselves both return a three-column matrix, consisting of the potential solution (xsol), the optimal solution in each projection plot (xopt), and the relative difference between the two (R = (xopt - xsol)/|xsol|).
An object of class summary.optproj inheriting from summary.optcheck, with elements:
xsolThe potential solution vector.
ysolThe value of the objective function at xsol.
maximizeLogical indicating whether the potential solution should maximize or minimize the objective function.
xoptA vector containing the argmax/argmin in each projection plot.
yoptA vector containing the max/min in each projection plot.
xdiffA two-column matrix containing the differences between xsol and xopt. The first column is the absolute difference D = xopt - xsol, the second is the relative difference R = D/|xsol|.
ydiffSame thing, but between ysol and yopt.
print.summary.optproj() for print method.
summary method for optrefit objects.summary method for optrefit objects.
## S3 method for class 'optrefit' summary(object, xnames, ...)## S3 method for class 'optrefit' summary(object, xnames, ...)
object |
An |
xnames |
Optional vector of names for the elements of the potential solution. |
... |
Further arguments to be passed to or from other methods. |
An object of class summary.optrefit inheriting from summary.optcheck, with elements:
xsolThe potential solution vector.
ysolThe value of the objective function at xsol.
maximizeLogical indicating whether the potential solution should maximize or minimize the objective function.
xoptA vector containing the argmax/argmin in each projection plot.
yoptThe scalar value of the max/min found by optim_refit.
xdiffA two-column matrix containing the differences between xsol and xopt. The first column is the absolute difference D = xopt - xsol, the second is the relative difference R = D/|xsol|.
ydiffA length-two vector containing the absolute and relative difference between ysol and yopt.
print.summary.optcheck() for print method.