find_extremum: Find Extremum of Fitted Lagrangian Multiplier Smoothing Spline

Description

Finds global extrema of a fitted lgspline model using deterministic or stochastic optimization strategies. Supports custom objective functions for advanced applications like Bayesian optimization acquisition functions.

Usage

find_extremum(
  object,
  vars = NULL,
  quick_heuristic = TRUE,
  initial = NULL,
  B_predict = NULL,
  minimize = FALSE,
  stochastic = FALSE,
  stochastic_draw = function(mu, sigma, ...) {
     N <- length(mu)
     rnorm(N, mu,
    sigma)
 },
  sigmasq_predict = object$sigmasq_tilde,
  custom_objective_function = NULL,
  custom_objective_derivative = NULL,
  ...
)

Value

A list containing the following components:

t: Numeric vector of input values at the extremum.
y: Numeric value of the objective function at the extremum.

Arguments

object

A fitted lgspline model object containing partition information and fitted values

vars

Vector; A vector of numeric indices (or character variable names) of predictors to optimize for. If NULL (by default), all predictors will be optimized.

quick_heuristic

Logical; whether to search only the top-performing partition. When TRUE (default), optimizes within the best partition. When FALSE, initiates searches from all partition local maxima.

initial

Numeric vector; Optional initial values for optimization. Useful for fixing binary predictors or providing starting points. Default NULL

B_predict

Matrix; Optional custom coefficient list for prediction. Useful for posterior draws in Bayesian optimization. Default NULL

minimize

Logical; whether to find minimum instead of maximum. Default FALSE

stochastic

Logical; whether to add noise for stochastic optimization. Enables better exploration of the function space. Default FALSE

stochastic_draw

Function; Generates random noise/modifies predictions for stochastic optimization, analogous to posterior_predictive_draw. Takes three arguments:

mu: Vector of predicted values
sigma: Vector of standard deviations (square-root of sigmasq_tilde)
...: Additional arguments to pass through

Default rnorm(length(mu), mu, sigma)

sigmasq_predict

Numeric; Variance parameter for stochastic optimization. Controls the magnitude of random perturbations. Defaults to object$sigmasq_tilde

custom_objective_function

Function; Optional custom objective function for optimization. Takes arguments:

mu: Vector of predicted response values
sigma: Vector of standard deviations
y_best: Numeric; Best observed response value
...: Additional arguments passed through

Default NULL

custom_objective_derivative

Function; Optional gradient function for custom optimization objective. Takes arguments:

mu: Vector of predicted response values
sigma: Vector of standard deviations
y_best: Numeric; Best observed response value
d_mu: Gradient of fitted function (for chain-rule computations)
...: Additional arguments passed through

Default NULL

...

Additional arguments passed to internal optimization routines.

Details

This method finds extrema (maxima or minima) of the fitted function or composite functions of the fit. The optimization process can be customized through several approaches:

Partition-based search: Either focuses on the top-performing partition (quick_heuristic = TRUE) or searches across all partition local maxima
Stochastic optimization: Adds random noise during optimization for better exploration
Custom objectives: Supports user-defined objective functions and gradients for specialized optimization tasks like Bayesian optimization

Examples

Run this code


## Basic usage with simulated data
set.seed(1234)
t <- runif(1000, -10, 10)
y <- 2*sin(t) + -0.06*t^2 + rnorm(length(t))
model_fit <- lgspline(t, y)
plot(model_fit)

## Find global maximum and minimum
max_point <- find_extremum(model_fit)
min_point <- find_extremum(model_fit, minimize = TRUE)
abline(v = max_point$t, col = 'blue')  # Add maximum point
abline(v = min_point$t, col = 'red')   # Add minimum point

## Advanced usage: custom objective functions
# expected improvement acquisition function
ei_custom_objective_function = function(mu, sigma, y_best, ...) {
  d <- y_best - mu
  d * pnorm(d/sigma) + sigma * dnorm(d/sigma)
}
# derivative of ei
ei_custom_objective_derivative = function(mu, sigma, y_best, d_mu, ...) {
  d <- y_best - mu
  z <- d/sigma
  d_z <- -d_mu/sigma
  pnorm(z)*d_mu - d*dnorm(z)*d_z + sigma*z*dnorm(z)*d_z
}

## Single iteration of Bayesian optimization
post_draw <- generate_posterior(model_fit)
acq <- find_extremum(model_fit,
                     stochastic = TRUE,  # Enable stochastic exploration
                     B_predict = post_draw$post_draw_coefficients,
                     sigmasq_predict = post_draw$post_draw_sigmasq,
                     custom_objective_function = ei_custom_objective_function,
                     custom_objective_derivative = ei_custom_objective_derivative)
abline(v = acq$t, col = 'green')  # Add acquisition point