use runmat_builtins::{
BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinOutputMode,
BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
LogicalArray, StructValue, Tensor, Value,
};
use runmat_macros::runtime_builtin;
use crate::builtins::common::spec::{
BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
ReductionNaN, ResidencyPolicy, ShapeRequirements,
};
use crate::builtins::math::optim::common::{
call_function, initial_guess, lookup_option, option_f64, option_string, value_to_real_vector,
value_to_scalar, vector_to_value,
};
use crate::builtins::math::optim::type_resolvers::nonlinear_solve_type;
use crate::{build_runtime_error, BuiltinResult, RuntimeError};
const NAME: &str = "fminunc";
const ALGORITHM: &str = "quasi-newton bfgs";
const DEFAULT_TOL_X: f64 = 1.0e-6;
const DEFAULT_TOL_FUN: f64 = 1.0e-6;
const DEFAULT_MAX_ITER: usize = 400;
const DEFAULT_MAX_FUN_EVALS: usize = 40000;
const MAX_ITER_LIMIT: usize = 1_000_000;
const MAX_FUN_EVAL_LIMIT: usize = 10_000_000;
const MAX_DENSE_BFGS_VARIABLES: usize = 2048;
const C1: f64 = 1.0e-4;
const C2: f64 = 0.9;
const MAX_LINE_SEARCH_ITERS: usize = 24;
macro_rules! output_x {
() => {
BuiltinParamDescriptor {
name: "x",
ty: BuiltinParamType::NumericArray,
arity: BuiltinParamArity::Required,
default: None,
description: "Estimated local minimizer.",
}
};
}
macro_rules! output_fval {
() => {
BuiltinParamDescriptor {
name: "fval",
ty: BuiltinParamType::NumericScalar,
arity: BuiltinParamArity::Required,
default: None,
description: "Objective value at x.",
}
};
}
macro_rules! output_exitflag {
() => {
BuiltinParamDescriptor {
name: "exitflag",
ty: BuiltinParamType::NumericScalar,
arity: BuiltinParamArity::Required,
default: None,
description: "Convergence status code.",
}
};
}
macro_rules! output_output {
() => {
BuiltinParamDescriptor {
name: "output",
ty: BuiltinParamType::Any,
arity: BuiltinParamArity::Required,
default: None,
description: "Iteration/function-count metadata struct.",
}
};
}
macro_rules! output_grad {
() => {
BuiltinParamDescriptor {
name: "grad",
ty: BuiltinParamType::NumericArray,
arity: BuiltinParamArity::Required,
default: None,
description: "Gradient at x.",
}
};
}
macro_rules! output_hessian {
() => {
BuiltinParamDescriptor {
name: "hessian",
ty: BuiltinParamType::NumericArray,
arity: BuiltinParamArity::Required,
default: None,
description: "Final approximate Hessian matrix.",
}
};
}
const FMINUNC_OUTPUT_X: [BuiltinParamDescriptor; 1] = [output_x!()];
const FMINUNC_OUTPUT_X_FVAL: [BuiltinParamDescriptor; 2] = [output_x!(), output_fval!()];
const FMINUNC_OUTPUT_X_FVAL_EXITFLAG: [BuiltinParamDescriptor; 3] =
[output_x!(), output_fval!(), output_exitflag!()];
const FMINUNC_OUTPUT_X_FVAL_EXITFLAG_OUTPUT: [BuiltinParamDescriptor; 4] = [
output_x!(),
output_fval!(),
output_exitflag!(),
output_output!(),
];
const FMINUNC_OUTPUT_X_FVAL_EXITFLAG_OUTPUT_GRAD: [BuiltinParamDescriptor; 5] = [
output_x!(),
output_fval!(),
output_exitflag!(),
output_output!(),
output_grad!(),
];
const FMINUNC_OUTPUT_ALL: [BuiltinParamDescriptor; 6] = [
output_x!(),
output_fval!(),
output_exitflag!(),
output_output!(),
output_grad!(),
output_hessian!(),
];
macro_rules! input_fun {
() => {
BuiltinParamDescriptor {
name: "fun",
ty: BuiltinParamType::Any,
arity: BuiltinParamArity::Required,
default: None,
description: "Scalar objective callback.",
}
};
}
macro_rules! input_x0 {
() => {
BuiltinParamDescriptor {
name: "x0",
ty: BuiltinParamType::Any,
arity: BuiltinParamArity::Required,
default: None,
description: "Initial guess scalar/vector/array.",
}
};
}
const FMINUNC_INPUTS_CORE: [BuiltinParamDescriptor; 2] = [input_fun!(), input_x0!()];
const FMINUNC_INPUTS_WITH_OPTIONS: [BuiltinParamDescriptor; 3] = [
input_fun!(),
input_x0!(),
BuiltinParamDescriptor {
name: "options",
ty: BuiltinParamType::Any,
arity: BuiltinParamArity::Optional,
default: None,
description: "Options struct from optimoptions or optimset.",
},
];
const FMINUNC_SIGNATURES: [BuiltinSignatureDescriptor; 12] = [
BuiltinSignatureDescriptor {
label: "x = fminunc(fun, x0)",
inputs: &FMINUNC_INPUTS_CORE,
outputs: &FMINUNC_OUTPUT_X,
},
BuiltinSignatureDescriptor {
label: "x = fminunc(fun, x0, options)",
inputs: &FMINUNC_INPUTS_WITH_OPTIONS,
outputs: &FMINUNC_OUTPUT_X,
},
BuiltinSignatureDescriptor {
label: "[x, fval] = fminunc(fun, x0)",
inputs: &FMINUNC_INPUTS_CORE,
outputs: &FMINUNC_OUTPUT_X_FVAL,
},
BuiltinSignatureDescriptor {
label: "[x, fval] = fminunc(fun, x0, options)",
inputs: &FMINUNC_INPUTS_WITH_OPTIONS,
outputs: &FMINUNC_OUTPUT_X_FVAL,
},
BuiltinSignatureDescriptor {
label: "[x, fval, exitflag] = fminunc(fun, x0)",
inputs: &FMINUNC_INPUTS_CORE,
outputs: &FMINUNC_OUTPUT_X_FVAL_EXITFLAG,
},
BuiltinSignatureDescriptor {
label: "[x, fval, exitflag] = fminunc(fun, x0, options)",
inputs: &FMINUNC_INPUTS_WITH_OPTIONS,
outputs: &FMINUNC_OUTPUT_X_FVAL_EXITFLAG,
},
BuiltinSignatureDescriptor {
label: "[x, fval, exitflag, output] = fminunc(fun, x0)",
inputs: &FMINUNC_INPUTS_CORE,
outputs: &FMINUNC_OUTPUT_X_FVAL_EXITFLAG_OUTPUT,
},
BuiltinSignatureDescriptor {
label: "[x, fval, exitflag, output] = fminunc(fun, x0, options)",
inputs: &FMINUNC_INPUTS_WITH_OPTIONS,
outputs: &FMINUNC_OUTPUT_X_FVAL_EXITFLAG_OUTPUT,
},
BuiltinSignatureDescriptor {
label: "[x, fval, exitflag, output, grad] = fminunc(fun, x0)",
inputs: &FMINUNC_INPUTS_CORE,
outputs: &FMINUNC_OUTPUT_X_FVAL_EXITFLAG_OUTPUT_GRAD,
},
BuiltinSignatureDescriptor {
label: "[x, fval, exitflag, output, grad] = fminunc(fun, x0, options)",
inputs: &FMINUNC_INPUTS_WITH_OPTIONS,
outputs: &FMINUNC_OUTPUT_X_FVAL_EXITFLAG_OUTPUT_GRAD,
},
BuiltinSignatureDescriptor {
label: "[x, fval, exitflag, output, grad, hessian] = fminunc(fun, x0)",
inputs: &FMINUNC_INPUTS_CORE,
outputs: &FMINUNC_OUTPUT_ALL,
},
BuiltinSignatureDescriptor {
label: "[x, fval, exitflag, output, grad, hessian] = fminunc(fun, x0, options)",
inputs: &FMINUNC_INPUTS_WITH_OPTIONS,
outputs: &FMINUNC_OUTPUT_ALL,
},
];
const FMINUNC_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
code: "RM.FMINUNC.INVALID_ARGUMENT",
identifier: Some("RunMat:fminunc:InvalidArgument"),
when: "Argument grammar/options parsing is invalid.",
message: "fminunc: invalid argument",
};
const FMINUNC_ERROR_INVALID_INPUT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
code: "RM.FMINUNC.INVALID_INPUT",
identifier: Some("RunMat:fminunc:InvalidInput"),
when: "Initial guess/callback/iteration semantics are invalid.",
message: "fminunc: invalid input",
};
const FMINUNC_ERRORS: [BuiltinErrorDescriptor; 2] =
[FMINUNC_ERROR_INVALID_ARGUMENT, FMINUNC_ERROR_INVALID_INPUT];
pub const FMINUNC_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
signatures: &FMINUNC_SIGNATURES,
output_mode: BuiltinOutputMode::ByRequestedOutputCount,
completion_policy: BuiltinCompletionPolicy::Public,
errors: &FMINUNC_ERRORS,
};
fn fminunc_error_with_detail(
error: &'static BuiltinErrorDescriptor,
detail: impl AsRef<str>,
) -> RuntimeError {
let detail = detail.as_ref();
let message = if detail.starts_with("fminunc:") {
detail.to_string()
} else {
format!("{}: {detail}", error.message)
};
let mut builder = build_runtime_error(message).with_builtin(NAME);
if let Some(identifier) = error.identifier {
builder = builder.with_identifier(identifier);
}
builder.build()
}
fn fminunc_map_error(err: RuntimeError, fallback: &'static BuiltinErrorDescriptor) -> RuntimeError {
if err.identifier().is_some() {
err
} else {
fminunc_error_with_detail(fallback, err.message())
}
}
#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::math::optim::fminunc")]
pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
name: "fminunc",
op_kind: GpuOpKind::Custom("unconstrained-min"),
supported_precisions: &[],
broadcast: BroadcastSemantics::None,
provider_hooks: &[],
constant_strategy: ConstantStrategy::InlineLiteral,
residency: ResidencyPolicy::GatherImmediately,
nan_mode: ReductionNaN::Include,
two_pass_threshold: None,
workgroup_size: None,
accepts_nan_mode: false,
notes: "Host BFGS solver with strong-Wolfe line search. Callback computations may use GPU-aware builtins, but optimizer state is gathered on the CPU.",
};
#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::math::optim::fminunc")]
pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
name: "fminunc",
shape: ShapeRequirements::Any,
constant_strategy: ConstantStrategy::InlineLiteral,
elementwise: None,
reduction: None,
emits_nan: false,
notes: "Unconstrained minimization repeatedly invokes user callbacks and terminates fusion planning.",
};
#[runtime_builtin(
name = "fminunc",
category = "math/optim",
summary = "Find an unconstrained local minimum of a smooth scalar objective.",
keywords = "fminunc,unconstrained minimization,bfgs,quasi-newton,strong wolfe,optimization",
accel = "sink",
type_resolver(nonlinear_solve_type),
descriptor(crate::builtins::math::optim::fminunc::FMINUNC_DESCRIPTOR),
builtin_path = "crate::builtins::math::optim::fminunc"
)]
async fn fminunc_builtin(function: Value, x0: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
if rest.len() > 1 {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
"too many input arguments",
));
}
let options_struct = parse_options(rest.first())
.map_err(|err| fminunc_map_error(err, &FMINUNC_ERROR_INVALID_ARGUMENT))?;
let guess = initial_guess(NAME, x0)
.await
.map_err(|err| fminunc_map_error(err, &FMINUNC_ERROR_INVALID_INPUT))?;
let options = FminuncOptions::from_struct(options_struct.as_ref(), guess.values.len())
.map_err(|err| fminunc_map_error(err, &FMINUNC_ERROR_INVALID_ARGUMENT))?;
let outcome = minimize(
&function,
guess.values,
&guess.shape,
guess.scalar,
&options,
)
.await
.map_err(|err| fminunc_map_error(err, &FMINUNC_ERROR_INVALID_INPUT))?;
finalize(outcome, &guess.shape, guess.scalar, &options)
}
fn parse_options(value: Option<&Value>) -> BuiltinResult<Option<StructValue>> {
match value {
None => Ok(None),
Some(Value::Struct(options)) => Ok(Some(options.clone())),
Some(other) => Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
format!("options must be a struct, got {other:?}"),
)),
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum DisplayMode {
Off,
Iter,
Notify,
Final,
}
impl DisplayMode {
fn parse(text: &str) -> BuiltinResult<Self> {
match text.to_ascii_lowercase().as_str() {
"off" | "none" => Ok(Self::Off),
"iter" => Ok(Self::Iter),
"notify" => Ok(Self::Notify),
"final" => Ok(Self::Final),
other => Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
format!(
"option Display must be 'off', 'iter', 'notify', or 'final', got '{other}'"
),
)),
}
}
}
#[derive(Debug, Clone)]
struct FminuncOptions {
tol_x: f64,
tol_fun: f64,
max_iter: usize,
max_fun_evals: usize,
display: DisplayMode,
specify_objective_gradient: bool,
}
impl FminuncOptions {
fn from_struct(options: Option<&StructValue>, variables: usize) -> BuiltinResult<Self> {
let display = DisplayMode::parse(&option_string(options, "Display", "off")?)?;
let algorithm = option_string(options, "Algorithm", "quasi-newton")?;
if !matches!(algorithm.as_str(), "quasi-newton" | "bfgs") {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
"fminunc: only Algorithm='quasi-newton' is supported",
));
}
let tol_x = option_f64(NAME, options, "TolX", DEFAULT_TOL_X)?;
let tol_fun = option_f64(NAME, options, "TolFun", DEFAULT_TOL_FUN)?;
if tol_x <= 0.0 || tol_fun <= 0.0 {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
"options TolX and TolFun must be positive",
));
}
let max_iter =
bounded_option_usize(options, "MaxIter", DEFAULT_MAX_ITER, MAX_ITER_LIMIT)?.max(1);
let max_fun_evals = bounded_option_usize(
options,
"MaxFunEvals",
DEFAULT_MAX_FUN_EVALS,
MAX_FUN_EVAL_LIMIT,
)?
.max(1);
let specify_objective_gradient = option_bool(options, "SpecifyObjectiveGradient", false)?;
validate_problem_size(variables)?;
Ok(Self {
tol_x,
tol_fun,
max_iter,
max_fun_evals,
display,
specify_objective_gradient,
})
}
}
fn bounded_option_usize(
options: Option<&StructValue>,
field: &str,
default: usize,
maximum: usize,
) -> BuiltinResult<usize> {
let value = option_f64(NAME, options, field, default as f64)?;
if value < 0.0 {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
format!("option {field} must be non-negative"),
));
}
if value > maximum as f64 {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
format!("option {field} must be no greater than {maximum}"),
));
}
if value.fract() != 0.0 {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
format!("option {field} must be an integer"),
));
}
Ok(value.floor() as usize)
}
fn validate_problem_size(variables: usize) -> BuiltinResult<()> {
if variables > MAX_DENSE_BFGS_VARIABLES {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_INPUT,
format!(
"problem has {variables} variables; current dense BFGS implementation supports at most {MAX_DENSE_BFGS_VARIABLES}"
),
));
}
dense_len(variables)?;
Ok(())
}
fn option_bool(options: Option<&StructValue>, field: &str, default: bool) -> BuiltinResult<bool> {
let Some(options) = options else {
return Ok(default);
};
let Some(value) = lookup_option(options, field) else {
return Ok(default);
};
bool_value(field, value)
}
fn bool_value(field: &str, value: &Value) -> BuiltinResult<bool> {
match value {
Value::Bool(flag) => Ok(*flag),
Value::Num(n) => bool_from_number(field, *n),
Value::Int(i) => bool_from_number(field, i.to_f64()),
Value::LogicalArray(LogicalArray { data, .. }) if data.len() == 1 => Ok(data[0] != 0),
Value::Tensor(Tensor { data, .. }) if data.len() == 1 => bool_from_number(field, data[0]),
Value::String(s) => bool_from_text(field, s),
Value::StringArray(sa) if sa.data.len() == 1 => bool_from_text(field, &sa.data[0]),
Value::CharArray(chars) if chars.rows == 1 => {
let text: String = chars.data.iter().collect();
bool_from_text(field, &text)
}
other => Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
format!("option {field} must be logical, got {other:?}"),
)),
}
}
fn bool_from_number(field: &str, value: f64) -> BuiltinResult<bool> {
if !value.is_finite() {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
format!("option {field} must be finite"),
));
}
if value == 0.0 {
Ok(false)
} else if value == 1.0 {
Ok(true)
} else {
Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
format!("option {field} must be logical 0 or 1"),
))
}
}
fn bool_from_text(field: &str, value: &str) -> BuiltinResult<bool> {
match value.trim().to_ascii_lowercase().as_str() {
"on" | "true" | "yes" => Ok(true),
"off" | "false" | "no" => Ok(false),
other => Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_ARGUMENT,
format!("option {field} must be 'on' or 'off', got '{other}'"),
)),
}
}
#[derive(Debug, Clone)]
struct ObjectiveState {
x: Vec<f64>,
f: f64,
grad: Vec<f64>,
}
#[derive(Debug, Clone)]
struct MinimizeOutcome {
x: Vec<f64>,
fval: f64,
grad: Vec<f64>,
inverse_hessian: Vec<f64>,
iterations: usize,
func_count: usize,
exitflag: i32,
firstorderopt: f64,
stepsize: f64,
message: String,
}
struct ObjectiveEvaluator<'a> {
function: &'a Value,
shape: Vec<usize>,
scalar: bool,
specify_gradient: bool,
max_fun_evals: usize,
func_count: usize,
}
impl<'a> ObjectiveEvaluator<'a> {
fn new(
function: &'a Value,
shape: &[usize],
scalar: bool,
specify_gradient: bool,
max_fun_evals: usize,
) -> Self {
Self {
function,
shape: shape.to_vec(),
scalar,
specify_gradient,
max_fun_evals,
func_count: 0,
}
}
async fn evaluate(&mut self, x: &[f64]) -> BuiltinResult<ObjectiveState> {
self.ensure_state_budget(x.len())?;
if self.specify_gradient {
self.evaluate_with_gradient(x).await
} else {
let f = self.evaluate_value(x).await?;
let grad = self.forward_difference_gradient(x, f).await?;
Ok(ObjectiveState {
x: x.to_vec(),
f,
grad,
})
}
}
async fn evaluate_value(&mut self, x: &[f64]) -> BuiltinResult<f64> {
self.ensure_value_budget()?;
let arg = vector_to_value(NAME, x.to_vec(), &self.shape, self.scalar)?;
let value = call_function(self.function, vec![arg]).await?;
let value = crate::dispatcher::gather_if_needed_async(&value).await?;
self.func_count += 1;
value_to_scalar(NAME, value)
}
async fn evaluate_with_gradient(&mut self, x: &[f64]) -> BuiltinResult<ObjectiveState> {
let callback =
crate::canonicalize_callback_handle_for_semantic_resolution(self.function.clone());
let arg = vector_to_value(NAME, x.to_vec(), &self.shape, self.scalar)?;
let value = crate::call_feval_async_with_outputs(callback, &[arg], 2).await?;
let value = crate::dispatcher::gather_if_needed_async(&value).await?;
self.func_count += 1;
let Value::OutputList(outputs) = value else {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_INPUT,
"objective must return [f, g] when SpecifyObjectiveGradient is true",
));
};
if outputs.len() < 2 {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_INPUT,
"objective must return both objective value and gradient",
));
}
let f = value_to_scalar(NAME, outputs[0].clone())?;
let grad = value_to_real_vector(NAME, outputs[1].clone()).await?;
if grad.len() != x.len() {
return Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_INPUT,
format!(
"objective gradient length {} does not match x length {}",
grad.len(),
x.len()
),
));
}
Ok(ObjectiveState {
x: x.to_vec(),
f,
grad,
})
}
async fn forward_difference_gradient(&mut self, x: &[f64], f0: f64) -> BuiltinResult<Vec<f64>> {
let mut grad = vec![0.0; x.len()];
for i in 0..x.len() {
let step = finite_difference_step(x[i]);
let mut trial = x.to_vec();
trial[i] += step;
let f_step = self.evaluate_value(&trial).await?;
grad[i] = (f_step - f0) / step;
}
Ok(grad)
}
fn state_evaluation_cost(&self, variables: usize) -> Option<usize> {
if self.specify_gradient {
Some(1)
} else {
variables.checked_add(1)
}
}
fn remaining_func_evals(&self) -> usize {
self.max_fun_evals.saturating_sub(self.func_count)
}
fn can_evaluate_state(&self, variables: usize) -> bool {
self.state_evaluation_cost(variables)
.is_some_and(|cost| cost <= self.remaining_func_evals())
}
fn can_evaluate_value(&self) -> bool {
self.remaining_func_evals() >= 1
}
fn ensure_state_budget(&self, variables: usize) -> BuiltinResult<()> {
if self.can_evaluate_state(variables) {
return Ok(());
}
Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_INPUT,
format!(
"MaxFunEvals exceeded before completing objective/gradient evaluation (FuncCount {}, MaxFunEvals {})",
self.func_count, self.max_fun_evals
),
))
}
fn ensure_value_budget(&self) -> BuiltinResult<()> {
if self.can_evaluate_value() {
return Ok(());
}
Err(fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_INPUT,
format!(
"MaxFunEvals exceeded before objective evaluation (FuncCount {}, MaxFunEvals {})",
self.func_count, self.max_fun_evals
),
))
}
}
fn finite_difference_step(x: f64) -> f64 {
f64::EPSILON.sqrt() * x.abs().max(1.0)
}
async fn minimize(
function: &Value,
x0: Vec<f64>,
shape: &[usize],
scalar: bool,
options: &FminuncOptions,
) -> BuiltinResult<MinimizeOutcome> {
let n = x0.len();
validate_problem_size(n)?;
let mut evaluator = ObjectiveEvaluator::new(
function,
shape,
scalar,
options.specify_objective_gradient,
options.max_fun_evals,
);
let mut state = evaluator.evaluate(&x0).await?;
let mut inverse_hessian = identity(n)?;
let mut iterations = 0usize;
let mut stepsize = 0.0f64;
let mut exitflag = 0i32;
let mut line_search_failed = false;
if norm_inf(&state.grad) <= options.tol_fun {
exitflag = 1;
}
while exitflag == 0
&& iterations < options.max_iter
&& evaluator.func_count < options.max_fun_evals
{
let mut direction = mat_vec(&inverse_hessian, &state.grad, n);
for value in &mut direction {
*value = -*value;
}
if dot(&direction, &state.grad) >= -1.0e-14 || !all_finite(&direction) {
inverse_hessian = identity(n)?;
direction = state.grad.iter().map(|value| -value).collect();
}
let next = match line_search(&mut evaluator, &state, &direction).await? {
LineSearchResult::Step(next) => next,
LineSearchResult::BudgetExhausted => {
exitflag = 0;
break;
}
LineSearchResult::Failed => {
line_search_failed = true;
exitflag = -3;
break;
}
};
let s = subtract(&next.x, &state.x);
stepsize = norm_inf(&s);
let y = subtract(&next.grad, &state.grad);
update_inverse_hessian(&mut inverse_hessian, &s, &y, n);
let f_delta = (state.f - next.f).abs();
let x_scale = 1.0 + norm_inf(&state.x);
let f_scale = 1.0 + state.f.abs();
state = next;
iterations += 1;
if norm_inf(&state.grad) <= options.tol_fun {
exitflag = 1;
} else if stepsize <= options.tol_x * x_scale {
exitflag = 2;
} else if f_delta <= options.tol_fun * f_scale {
exitflag = 3;
}
}
let firstorderopt = norm_inf(&state.grad);
let message = build_message(
exitflag,
iterations,
evaluator.func_count,
line_search_failed,
);
emit_summary(
&state,
exitflag,
&message,
options,
iterations,
evaluator.func_count,
);
Ok(MinimizeOutcome {
x: state.x,
fval: state.f,
grad: state.grad,
inverse_hessian,
iterations,
func_count: evaluator.func_count,
exitflag,
firstorderopt,
stepsize,
message,
})
}
async fn line_search(
evaluator: &mut ObjectiveEvaluator<'_>,
state: &ObjectiveState,
direction: &[f64],
) -> BuiltinResult<LineSearchResult> {
let dphi0 = dot(&state.grad, direction);
if dphi0 >= 0.0 || !dphi0.is_finite() {
return Ok(LineSearchResult::Failed);
}
let mut alpha_prev = 0.0;
let mut f_prev = state.f;
let mut alpha = 1.0;
for iter in 0..MAX_LINE_SEARCH_ITERS {
if !evaluator.can_evaluate_state(state.x.len()) {
return Ok(LineSearchResult::BudgetExhausted);
}
let trial_x = add_scaled(&state.x, direction, alpha);
let trial = evaluator.evaluate(&trial_x).await?;
if trial.f > state.f + C1 * alpha * dphi0 || (iter > 0 && trial.f >= f_prev) {
return zoom(evaluator, state, direction, alpha_prev, alpha).await;
}
let dphi = dot(&trial.grad, direction);
if dphi.abs() <= -C2 * dphi0 {
return Ok(LineSearchResult::Step(trial));
}
if dphi >= 0.0 {
return zoom(evaluator, state, direction, alpha, alpha_prev).await;
}
alpha_prev = alpha;
f_prev = trial.f;
alpha = (alpha * 2.0).min(64.0);
}
Ok(LineSearchResult::Failed)
}
enum LineSearchResult {
Step(ObjectiveState),
BudgetExhausted,
Failed,
}
async fn zoom(
evaluator: &mut ObjectiveEvaluator<'_>,
state: &ObjectiveState,
direction: &[f64],
mut lo: f64,
mut hi: f64,
) -> BuiltinResult<LineSearchResult> {
let dphi0 = dot(&state.grad, direction);
let mut f_lo = if lo == 0.0 {
state.f
} else {
if !evaluator.can_evaluate_value() {
return Ok(LineSearchResult::BudgetExhausted);
}
let x_lo = add_scaled(&state.x, direction, lo);
evaluator.evaluate_value(&x_lo).await?
};
for _ in 0..MAX_LINE_SEARCH_ITERS {
if !evaluator.can_evaluate_state(state.x.len()) {
return Ok(LineSearchResult::BudgetExhausted);
}
let alpha = 0.5 * (lo + hi);
let trial_x = add_scaled(&state.x, direction, alpha);
let trial = evaluator.evaluate(&trial_x).await?;
if trial.f > state.f + C1 * alpha * dphi0 || trial.f >= f_lo {
hi = alpha;
} else {
let dphi = dot(&trial.grad, direction);
if dphi.abs() <= -C2 * dphi0 {
return Ok(LineSearchResult::Step(trial));
}
if dphi * (hi - lo) >= 0.0 {
hi = lo;
}
lo = alpha;
f_lo = trial.f;
}
if (hi - lo).abs() <= 1.0e-12 * (1.0 + lo.abs() + hi.abs()) {
if trial.f <= state.f + C1 * alpha * dphi0 {
return Ok(LineSearchResult::Step(trial));
}
return Ok(LineSearchResult::Failed);
}
}
Ok(LineSearchResult::Failed)
}
fn update_inverse_hessian(h: &mut [f64], s: &[f64], y: &[f64], n: usize) {
let ys = dot(y, s);
if ys <= 1.0e-14 || !ys.is_finite() {
return;
}
let rho = 1.0 / ys;
let hy = mat_vec(h, y, n);
let yhy = dot(y, &hy);
let coeff = (1.0 + yhy * rho) * rho;
for row in 0..n {
for col in 0..n {
h[row * n + col] +=
coeff * s[row] * s[col] - rho * (s[row] * hy[col] + hy[row] * s[col]);
}
}
}
fn finalize(
outcome: MinimizeOutcome,
shape: &[usize],
scalar: bool,
_options: &FminuncOptions,
) -> BuiltinResult<Value> {
let x = vector_to_value(NAME, outcome.x.clone(), shape, scalar)?;
let fval = Value::Num(outcome.fval);
let exitflag = Value::Num(outcome.exitflag as f64);
let output = Value::Struct(build_output_struct(&outcome));
let grad = vector_to_value(NAME, outcome.grad.clone(), shape, scalar)?;
let hessian = hessian_value(&outcome.inverse_hessian, outcome.grad.len(), scalar)?;
match crate::output_count::current_output_count() {
None => Ok(x),
Some(0) => Ok(Value::OutputList(Vec::new())),
Some(1) => Ok(crate::output_count::output_list_with_padding(1, vec![x])),
Some(2) => Ok(crate::output_count::output_list_with_padding(
2,
vec![x, fval],
)),
Some(3) => Ok(crate::output_count::output_list_with_padding(
3,
vec![x, fval, exitflag],
)),
Some(4) => Ok(crate::output_count::output_list_with_padding(
4,
vec![x, fval, exitflag, output],
)),
Some(5) => Ok(crate::output_count::output_list_with_padding(
5,
vec![x, fval, exitflag, output, grad],
)),
Some(n) if n >= 6 => Ok(crate::output_count::output_list_with_padding(
n,
vec![x, fval, exitflag, output, grad, hessian],
)),
Some(_) => unreachable!("output count cases above are exhaustive"),
}
}
fn hessian_value(inverse_hessian: &[f64], n: usize, scalar: bool) -> BuiltinResult<Value> {
let hessian = match invert_matrix(inverse_hessian, n) {
Some(hessian) => hessian,
None => identity(n)?,
};
if scalar {
Ok(Value::Num(hessian[0]))
} else {
Tensor::new(hessian, vec![n, n])
.map(Value::Tensor)
.map_err(|e| fminunc_error_with_detail(&FMINUNC_ERROR_INVALID_INPUT, e))
}
}
fn build_output_struct(outcome: &MinimizeOutcome) -> StructValue {
let mut fields = StructValue::new();
fields.insert("iterations", Value::Num(outcome.iterations as f64));
fields.insert("funcCount", Value::Num(outcome.func_count as f64));
fields.insert("algorithm", Value::from(ALGORITHM));
fields.insert("firstorderopt", Value::Num(outcome.firstorderopt));
fields.insert("stepsize", Value::Num(outcome.stepsize));
fields.insert("message", Value::from(outcome.message.clone()));
fields
}
fn build_message(
exitflag: i32,
iterations: usize,
func_count: usize,
line_search_failed: bool,
) -> String {
match exitflag {
1 => format!(
"Optimization terminated: first-order optimality is below OPTIONS.TolFun. Iterations: {iterations}, FuncCount: {func_count}."
),
2 => format!(
"Optimization terminated: step size is below OPTIONS.TolX. Iterations: {iterations}, FuncCount: {func_count}."
),
3 => format!(
"Optimization terminated: objective change is below OPTIONS.TolFun. Iterations: {iterations}, FuncCount: {func_count}."
),
-3 if line_search_failed => format!(
"Exiting: line search could not find a point satisfying the strong Wolfe conditions. Iterations: {iterations}, FuncCount: {func_count}."
),
_ => format!(
"Exiting: Maximum number of function evaluations or iterations has been exceeded. Iterations: {iterations}, FuncCount: {func_count}."
),
}
}
fn emit_summary(
state: &ObjectiveState,
exitflag: i32,
message: &str,
options: &FminuncOptions,
iterations: usize,
func_count: usize,
) {
let should_emit = match options.display {
DisplayMode::Off => false,
DisplayMode::Final | DisplayMode::Iter => true,
DisplayMode::Notify => exitflag <= 0,
};
if !should_emit {
return;
}
crate::console::record_console_line(
crate::console::ConsoleStream::Stdout,
format!(
"fminunc: fval = {fval:.6e}, firstorderopt = {opt:.6e}, exitflag = {exitflag}. {message}",
fval = state.f,
opt = norm_inf(&state.grad),
),
);
if matches!(options.display, DisplayMode::Iter) {
crate::console::record_console_line(
crate::console::ConsoleStream::Stdout,
format!("fminunc: iterations = {iterations}, funcCount = {func_count}"),
);
}
}
fn dense_len(n: usize) -> BuiltinResult<usize> {
n.checked_mul(n).ok_or_else(|| {
fminunc_error_with_detail(
&FMINUNC_ERROR_INVALID_INPUT,
"dense BFGS matrix size overflows usize",
)
})
}
fn identity(n: usize) -> BuiltinResult<Vec<f64>> {
let len = dense_len(n)?;
let mut out = vec![0.0; len];
for i in 0..n {
out[i * n + i] = 1.0;
}
Ok(out)
}
fn invert_matrix(matrix: &[f64], n: usize) -> Option<Vec<f64>> {
let mut a = matrix.to_vec();
let mut inv = identity(n).ok()?;
for col in 0..n {
let mut pivot = col;
let mut pivot_abs = a[col * n + col].abs();
for row in (col + 1)..n {
let candidate = a[row * n + col].abs();
if candidate > pivot_abs {
pivot = row;
pivot_abs = candidate;
}
}
if pivot_abs <= 1.0e-14 || !pivot_abs.is_finite() {
return None;
}
if pivot != col {
for j in 0..n {
a.swap(col * n + j, pivot * n + j);
inv.swap(col * n + j, pivot * n + j);
}
}
let diag = a[col * n + col];
for j in 0..n {
a[col * n + j] /= diag;
inv[col * n + j] /= diag;
}
for row in 0..n {
if row == col {
continue;
}
let factor = a[row * n + col];
if factor == 0.0 {
continue;
}
for j in 0..n {
a[row * n + j] -= factor * a[col * n + j];
inv[row * n + j] -= factor * inv[col * n + j];
}
}
}
Some(inv)
}
fn mat_vec(matrix: &[f64], vector: &[f64], n: usize) -> Vec<f64> {
let mut out = vec![0.0; n];
for row in 0..n {
for col in 0..n {
out[row] += matrix[row * n + col] * vector[col];
}
}
out
}
fn dot(a: &[f64], b: &[f64]) -> f64 {
a.iter().zip(b.iter()).map(|(x, y)| x * y).sum()
}
fn norm_inf(values: &[f64]) -> f64 {
values.iter().fold(0.0, |acc, value| acc.max(value.abs()))
}
fn all_finite(values: &[f64]) -> bool {
values.iter().all(|value| value.is_finite())
}
fn subtract(a: &[f64], b: &[f64]) -> Vec<f64> {
a.iter().zip(b.iter()).map(|(x, y)| x - y).collect()
}
fn add_scaled(x: &[f64], direction: &[f64], alpha: f64) -> Vec<f64> {
x.iter()
.zip(direction.iter())
.map(|(xi, di)| xi + alpha * di)
.collect()
}
#[cfg(test)]
mod tests {
use super::*;
use futures::executor::block_on;
use runmat_builtins::IntValue;
use std::sync::Arc;
fn bound(id: usize) -> Value {
Value::BoundFunctionHandle {
name: format!("objective_{id}"),
function: id,
}
}
fn tensor(data: Vec<f64>) -> Value {
Value::Tensor(Tensor::new(data, vec![3, 1]).unwrap())
}
fn vector_from_value(value: &Value) -> Vec<f64> {
match value {
Value::Tensor(tensor) => tensor.data.clone(),
Value::Num(n) => vec![*n],
other => panic!("expected numeric value, got {other:?}"),
}
}
fn assert_close_vec(actual: &[f64], expected: &[f64], tol: f64) {
assert_eq!(actual.len(), expected.len());
for (idx, (a, e)) in actual.iter().zip(expected.iter()).enumerate() {
assert!((a - e).abs() <= tol, "at {idx}: expected {e}, got {a}");
}
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_quadratic_reaches_vector_minimum() {
let _guard = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
|_function, args, requested_outputs| {
assert_eq!(requested_outputs, 1);
let x = vector_from_value(&args[0]);
Box::pin(async move {
let target = [1.0, 2.0, 3.0];
let f = x
.iter()
.zip(target.iter())
.map(|(xi, ti)| (xi - ti).powi(2))
.sum::<f64>();
Ok(Value::Num(f))
})
},
)));
let result = block_on(fminunc_builtin(
bound(101),
tensor(vec![0.0, 0.0, 0.0]),
Vec::new(),
))
.expect("fminunc");
assert_close_vec(&vector_from_value(&result), &[1.0, 2.0, 3.0], 1.0e-4);
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_rosenbrock_converges_from_standard_start() {
let _guard = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
|_function, args, requested_outputs| {
assert_eq!(requested_outputs, 1);
let x = vector_from_value(&args[0]);
Box::pin(async move {
let a = x[0];
let b = x[1];
Ok(Value::Num(100.0 * (b - a * a).powi(2) + (1.0 - a).powi(2)))
})
},
)));
let x0 = Value::Tensor(Tensor::new(vec![-1.2, 1.0], vec![2, 1]).unwrap());
let result = block_on(fminunc_builtin(bound(102), x0, Vec::new())).expect("fminunc");
assert_close_vec(&vector_from_value(&result), &[1.0, 1.0], 2.0e-3);
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_handles_high_dimensional_smooth_objective() {
let _guard = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
|_function, args, requested_outputs| {
assert_eq!(requested_outputs, 1);
let x = vector_from_value(&args[0]);
Box::pin(async move {
let f = x
.iter()
.enumerate()
.map(|(idx, xi)| {
let target = (idx + 1) as f64;
(xi - target).powi(2)
})
.sum::<f64>();
Ok(Value::Num(f))
})
},
)));
let x0 = Value::Tensor(Tensor::new(vec![0.0; 8], vec![8, 1]).unwrap());
let result = block_on(fminunc_builtin(bound(103), x0, Vec::new())).expect("fminunc");
assert_close_vec(
&vector_from_value(&result),
&[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0],
1.0e-4,
);
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_uses_objective_gradient_and_returns_six_outputs() {
let _guard = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
|_function, args, requested_outputs| {
assert_eq!(requested_outputs, 2);
let x = vector_from_value(&args[0]);
Box::pin(async move {
let target = [1.0, 2.0, 3.0];
let grad = x
.iter()
.zip(target.iter())
.map(|(xi, ti)| 2.0 * (xi - ti))
.collect::<Vec<_>>();
let f = x
.iter()
.zip(target.iter())
.map(|(xi, ti)| (xi - ti).powi(2))
.sum::<f64>();
Ok(Value::OutputList(vec![
Value::Num(f),
Value::Tensor(Tensor::new(grad, vec![3, 1]).unwrap()),
]))
})
},
)));
let mut opts = StructValue::new();
opts.insert("SpecifyObjectiveGradient", Value::Bool(true));
opts.insert("TolX", Value::Num(1.0e-10));
opts.insert("TolFun", Value::Num(1.0e-10));
let _outputs = crate::output_count::push_output_count(Some(6));
let result = block_on(fminunc_builtin(
bound(104),
tensor(vec![0.0, 0.0, 0.0]),
vec![Value::Struct(opts)],
))
.expect("fminunc");
let Value::OutputList(outputs) = result else {
panic!("expected output list");
};
assert_eq!(outputs.len(), 6);
assert_close_vec(&vector_from_value(&outputs[0]), &[1.0, 2.0, 3.0], 1.0e-7);
assert!(matches!(outputs[1], Value::Num(f) if f.abs() < 1.0e-14));
assert!(matches!(outputs[2], Value::Num(flag) if flag > 0.0));
match &outputs[3] {
Value::Struct(output) => {
assert!(matches!(
output.fields.get("iterations"),
Some(Value::Num(_))
));
assert!(matches!(
output.fields.get("funcCount"),
Some(Value::Num(_))
));
assert!(
matches!(output.fields.get("firstorderopt"), Some(Value::Num(v)) if *v < 1.0e-8)
);
assert!(
matches!(output.fields.get("algorithm"), Some(Value::String(text)) if text.contains("bfgs"))
);
}
other => panic!("unexpected output struct {other:?}"),
}
assert_close_vec(&vector_from_value(&outputs[4]), &[0.0, 0.0, 0.0], 1.0e-7);
match &outputs[5] {
Value::Tensor(hessian) => {
assert_eq!(hessian.shape, vec![3, 3]);
assert!(hessian.data.iter().all(|value| value.is_finite()));
assert!(hessian.data[0] > 0.0);
assert!(hessian.data[4] > 0.0);
assert!(hessian.data[8] > 0.0);
}
other => panic!("unexpected hessian {other:?}"),
}
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_rejects_invalid_gradient_length() {
let _guard = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
|_function, _args, requested_outputs| {
assert_eq!(requested_outputs, 2);
Box::pin(async {
Ok(Value::OutputList(vec![
Value::Num(1.0),
Value::Tensor(Tensor::new(vec![1.0, 2.0], vec![2, 1]).unwrap()),
]))
})
},
)));
let mut opts = StructValue::new();
opts.insert("SpecifyObjectiveGradient", Value::Bool(true));
let err = block_on(fminunc_builtin(
bound(105),
tensor(vec![0.0, 0.0, 0.0]),
vec![Value::Struct(opts)],
))
.unwrap_err();
assert!(err.message().contains("gradient length"));
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_scalar_input_returns_scalar() {
let _guard = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
|_function, args, requested_outputs| {
assert_eq!(requested_outputs, 1);
let x = match &args[0] {
Value::Num(value) => *value,
other => panic!("expected scalar, got {other:?}"),
};
Box::pin(async move { Ok(Value::Num((x - 4.0).powi(2))) })
},
)));
let result =
block_on(fminunc_builtin(bound(106), Value::Num(0.0), Vec::new())).expect("fminunc");
assert!(matches!(result, Value::Num(x) if (x - 4.0).abs() < 1.0e-4));
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_rejects_non_struct_options() {
let err = block_on(fminunc_builtin(
bound(107),
Value::Num(0.0),
vec![Value::Int(IntValue::I32(1))],
))
.unwrap_err();
assert_eq!(err.identifier(), Some("RunMat:fminunc:InvalidArgument"));
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_rejects_budget_too_small_for_finite_difference_gradient() {
let _guard = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
|_function, _args, _requested_outputs| {
Box::pin(async { panic!("budget check should happen before invoking objective") })
},
)));
let mut opts = StructValue::new();
opts.insert("MaxFunEvals", Value::Num(3.0));
let err = block_on(fminunc_builtin(
bound(108),
tensor(vec![0.0, 0.0, 0.0]),
vec![Value::Struct(opts)],
))
.unwrap_err();
assert!(err.message().contains("MaxFunEvals"));
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_rejects_problem_too_large_for_dense_bfgs() {
let variables = MAX_DENSE_BFGS_VARIABLES + 1;
let x0 = Value::Tensor(Tensor::new(vec![0.0; variables], vec![variables, 1]).unwrap());
let err = block_on(fminunc_builtin(bound(109), x0, Vec::new())).unwrap_err();
assert!(err.message().contains("dense BFGS"));
}
#[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
#[test]
fn fminunc_rejects_huge_direct_iteration_limits() {
let mut opts = StructValue::new();
opts.insert("MaxFunEvals", Value::Num((MAX_FUN_EVAL_LIMIT as f64) + 1.0));
let err = block_on(fminunc_builtin(
bound(110),
Value::Num(0.0),
vec![Value::Struct(opts)],
))
.unwrap_err();
assert!(err.message().contains("MaxFunEvals"));
}
}