optimize

package module
v0.2.0 Latest Latest
Warning

This package is not in the latest version of its module.

Go to latest
Published: Oct 17, 2026 License: BSD-3-Clause Imports: 9 Imported by: 0

README

go-numerics/optimize

optimize

Docs pkg.go.dev License

Numerical optimisation and root finding in pure Go (CGO=0): minimisation, scalar and system roots, nonlinear least squares, curve fitting with parameter covariance, and linear programming.

🌐 Website · 📚 Documentation

The package covers the parts of scipy.optimize that scientific notebooks use most. Each function documents the SciPy call it corresponds to. The test suite replays SciPy 1.18.1 on standard problems, and it checks solutions, objective values, status codes, messages, and iteration and evaluation counts.

go get github.com/go-numerics/optimize

It works with go-ndarray: curve-fit data, covariances and LP constraint matrices are *ndarray.Array, and results have …Array() accessors.

What is there

Go SciPy Algorithm
Minimize(p, x0, NelderMead, opts) minimize(method='Nelder-Mead') downhill simplex, adaptive option, bounds
Minimize(…, Powell, …) minimize(method='Powell') conjugate directions, bounds
Minimize(…, CG, …) minimize(method='CG') Polak–Ribière+, strong Wolfe
Minimize(…, BFGS, …) minimize(method='BFGS') inverse-Hessian BFGS, O(n²) update
Minimize(…, NewtonCG, …) minimize(method='Newton-CG') truncated Newton (Hess, HessP or FD)
Minimize(…, LBFGSB, …) minimize(method='L-BFGS-B') L-BFGS-B 3.0 with bounds
MinimizeScalar, Fminbound, Bracket minimize_scalar (brent, golden, bounded), fminbound, bracket Brent, golden section, FMIN
Brentq, Brenth, Ridder, Bisect, Newton, RootScalar brentq, brenth, ridder, bisect, newton (Newton/secant/Halley), root_scalar
Root, Fsolve root(method='hybr'), fsolve MINPACK hybrd/hybrj (Powell hybrid)
LeastSquares(…, TRF, …) least_squares(method='trf') trust-region reflective, bounds, robust losses, x_scale
LeastSquares(…, LM, …), Leastsq least_squares(method='lm'), leastsq (with cov_x) MINPACK lmder/lmdif
CurveFit, CurveFitModel curve_fit (sigma 1-D/2-D, absolute_sigma, bounds, losses) LM or TRF, covariance as SciPy
Linprog linprog revised simplex (default) or interior point; duals
Brute brute grid then Nelder–Mead polish
CheckGrad, ApproxFprime check_grad, approx_fprime forward differences

Gradients are optional everywhere. Without one, finite differences use the step rules of SciPy's approx_derivative, kept inside the bounds. Non-convergence is reported in Success/Status/Message as SciPy does. Invalid input returns an error.

Examples

Each example below is a runnable, output-checked Go example in example_test.go, and its output was checked against SciPy.

// scipy.optimize.minimize(rosen, x0, jac=rosen_der, method='BFGS')
r, err := optimize.Minimize(optimize.Problem{Func: rosen, Grad: rosenDer},
	[]float64{-1.2, 1}, optimize.BFGS, nil)
// r.X = [1 1], r.Nit = 32 (SciPy: 32), "Optimization terminated successfully."

// scipy.optimize.minimize(f, x0, method='L-BFGS-B', bounds=[(-2, 0.5), (None, None)])
r, _ = optimize.Minimize(optimize.Problem{Func: f}, []float64{0, 0}, optimize.LBFGSB,
	&optimize.Options{Bounds: []optimize.Bound{{Lo: -2, Hi: 0.5}, optimize.Unbounded()}})

// popt, pcov = scipy.optimize.curve_fit(lambda x, a, b, c: a*np.exp(-b*x) + c, x, y)
fit, err := optimize.CurveFit(func(x float64, p []float64) float64 {
	return p[0]*math.Exp(-p[1]*x) + p[2]
}, xdata, ydata, []float64{1, 1, 1}, nil)
// fit.Popt, fit.PcovArray(), fit.Perr()  ->  a = 2.500 ± 0.004, b = 1.300 ± 0.005, c = 0.500 ± 0.002

// scipy.optimize.brentq(lambda x: x**3 - 2*x - 5, 2, 3)
root, _ := optimize.Brentq(func(x float64) float64 { return x*x*x - 2*x - 5 }, 2, 3, nil)
// root.Root = 2.094551481542…, 7 iterations (as SciPy)

// scipy.optimize.root(fun, [0, 0])        (fun writes f(x) into its second argument)
sol, _ := optimize.Root(optimize.SystemProblem{Func: fun}, []float64{0, 0}, nil)

// scipy.optimize.least_squares(fun, x0, bounds=([0, 0], [2, 0.5]), loss='soft_l1', f_scale=0.1)
ls, _ := optimize.LeastSquares(optimize.LSProblem{Func: resid, M: m}, []float64{1, 0.1}, optimize.TRF,
	&optimize.LSOptions{Bounds: []optimize.Bound{{0, 2}, {0, 0.5}}, Loss: optimize.LossSoftL1, FScale: 0.1})

// scipy.optimize.linprog(c, A_ub=A, b_ub=b, bounds=[(None, None), (-3, None)])
lp, _ := optimize.Linprog([]float64{-1, 4}, optimize.LPConstraints{AUb: A, BUb: b,
	Bounds: []optimize.Bound{optimize.Unbounded(), {Lo: -3, Hi: math.Inf(1)}}}, nil)
// lp.X = [10 -3], lp.Fun = -22, lp.IneqLin.Marginals (HiGHS sign convention)

Callbacks receive the point first and write their outputs into the slices that follow: Func(x), Grad(x, g), FuncGrad(x, g) (SciPy's jac=True), Hess(x, h), HessP(x, p, hp), residuals Func(x, f), Jacobians Jac(x, J). Matrices are row-major []float64.

How close to SciPy

All numbers below are against SciPy 1.18.1 fixtures generated by testdata/gen/*.py (run with uv run --with scipy --with numpy python …) and embedded in the tests. They are the same on arm64 and amd64, except where noted.

Area Fixtures Result
Minimize (NM, Powell, CG, BFGS, Newton-CG, L-BFGS-B) 79 runs on Rosenbrock n=2/10/100, Beale, Himmelblau, Powell singular, Wood, with and without gradients, bounds, iteration limits every converged SciPy run is matched with an objective at least as good, within the method's own tolerance; on the arithmetic-only problems, 36/38 runs take SciPy's exact path (identical nit, nfev, njev); the 2 exceptions are explained in docs/design.md
L-BFGS-B 44 runs, including 8 problems that reach rare branches (Cholesky failures, the 3.0 backtrack, Cauchy breakpoints) identical nit/nfev/njev/status/message on every analytic-gradient run except Rosenbrock n=100 (513 vs 514 iterations, SciPy's Accelerate BLAS rounding); bit-identical to SciPy's C code built with sequential BLAS
Root/Fsolve 23 runs (MGH problems with and without Jacobian, failure and option cases) 23/23 identical counts, status and message; x bit-identical
LeastSquares(LM) / Leastsq 39 / 44 runs on the MGH least-squares set 36/39 and 41/44 identical counts; the mismatches are on problems that call exp/sin/cos, where Go's libm and the C library differ in the last ulp
LeastSquares(TRF) 89 runs: MGH problems, bounds, all robust losses, x_scale, 2-/3-point differences 86/89 identical nfev/njev/status (84/89 on amd64); every cost and active mask matches
CurveFit 11 fits (exponentials, Gaussians, sum of two Gaussians, a 2-variable plane, 1-D and 2-D sigma, bounds, absolute_sigma, m = n) popt within 1e-6, pcov within 1e-5 relative, identical nfev on every fit
Scalar roots 71 runs on 11 hard functions (multiple roots, flat, steep, kinks, poles) × 7 methods 68/71 on SciPy's exact path (same root, iterations, calls); the rest differ only through cos/exp in the last ulp
MinimizeScalar, Bracket 36 + 6 runs identical x/nit/nfev/message on arithmetic-only functions
Linprog 46 problems vs HiGHS: Klee–Minty, cycling examples, degenerate, free and upper-bounded variables, redundant rows, infeasible, unbounded, random up to 60×40 46/46 statuses; objective error ≤ 1.8e-15 (simplex) and ≤ 2e-8 (interior point); solutions and duals inside SciPy's optimal sets
Brute 4 grids (the SciPy documentation example, Himmelblau, Rosenbrock 3-D) x and f to 1e-9

How such close agreement is possible:

  • The ports follow their references statement by statement.
  • Every product that feeds an addition is rounded explicitly, because Go would otherwise fuse it into an FMA on arm64 and other architectures.
  • The tests use FMA-free oracles.

See docs/design.md.

Negative controls. One algorithmic step was broken at a time; each break made the suite fail, and each was then reverted:

Broken step Failing tests
A term dropped from the BFGS update 6
A sign flipped in Moré–Thuente's dcstep 57
The wrong Brent golden constant 4
The Levenberg correction dropped in lmpar 40
A sign flipped in dogleg 17
TRF reflection disabled, or its radius growth changed 20 / 48 of the 89 cases leave SciPy's path
Max instead of min in the simplex ratio test 83
The sign of the interior-point corrector flipped iteration parity drops from 43 to 2 of 46
L-BFGS-B's subspace backtrack skipped, its Cauchy f2 update dropped, or its breakpoint heap disabled named fixtures fail for each

Performance

Medians of 7 rounds, with Go and SciPy interleaved on a loaded host. docs/perf.md has the method and caveats.

Benchmark Go SciPy Speed-up
Rosenbrock n=2, BFGS 5.7 µs 2.58 ms 449×
Rosenbrock n=100, BFGS 14.8 ms 87.9 ms 5.9×
Rosenbrock n=2, L-BFGS-B 126 µs 1.47 ms 11.7×
Rosenbrock n=100, L-BFGS-B 6.36 ms 23.2 ms 3.7×
curve_fit, 10⁴ points, LM 4.13 ms 2.78 ms 0.67× (SciPy faster)
curve_fit, 10⁴ points, TRF 5.97 ms 4.61 ms 0.77× (SciPy faster)
brentq 0.27 µs 8.6 µs 32×

Most of SciPy's time on small problems is Python overhead around the callbacks. Even after subtracting SciPy's time spent in its own callbacks, the Go algorithms remain 2–300× faster, depending on the problem.

On large vectorised curve fits SciPy is faster today. Its model call is one NumPy expression, and its MINPACK/LAPACK kernels are compiled C. Our model evaluation is at parity; the gap is in the solver's workspace handling.

Not done

  • minimize: trust-constr, SLSQP, COBYLA/COBYQA, TNC, dogleg/trust-ncg and the other trust-region methods, and constraints other than bounds. For L-BFGS-B, the hess_inv operator is not provided.
  • root / least_squares: root methods other than hybr; the band/col_deriv options; dogbox; the lsmr trust-region solver and sparse Jacobians.
  • linprog: HiGHS-class performance, MILP, singleton-row presolve, crossover and sparse matrices.
  • Global optimisers: differential_evolution, basinhopping, shgo, dual_annealing.

Licence

BSD 3-Clause; see LICENSE.

Documentation

Overview

Package optimize provides numerical optimisation and root finding in pure Go: unconstrained and bound-constrained minimisation (Nelder–Mead, Powell, CG, BFGS, Newton-CG, L-BFGS-B), scalar minimisation (Brent, golden section, bounded), scalar roots (Brent, Ridder, bisection, Newton, secant, Halley), systems of equations (MINPACK's hybrid Powell method), nonlinear least squares (MINPACK Levenberg–Marquardt and trust-region reflective with bounds and robust losses), curve fitting with parameter covariance, and linear programming.

The API follows scipy.optimize closely enough that each function documents the SciPy call it corresponds to, and the test suite checks the results against SciPy fixtures; the implementation is independent Go.

Index

Examples

Constants

View Source
const (
	LPOptimal       = 0 // optimisation terminated successfully
	LPIterLimit     = 1 // iteration limit reached
	LPInfeasible    = 2 // problem appears to be infeasible
	LPUnbounded     = 3 // problem appears to be unbounded
	LPNumericalFail = 4 // numerical difficulties encountered
)

LP status codes, as scipy.optimize.linprog reports them in status.

View Source
const (
	FlagConverged = "converged"
	FlagConvErr   = "convergence error"
	FlagSignErr   = "sign error"
)

Root-finder flags, as SciPy's RootResults.flag strings.

Variables

View Source
var (
	// ErrDimension reports slices whose lengths disagree.
	ErrDimension = errors.New("optimize: dimension mismatch")
	// ErrMethod reports an unknown or unsupported method.
	ErrMethod = errors.New("optimize: unknown method")
	// ErrBounds reports inconsistent bounds or a starting point outside them.
	ErrBounds = errors.New("optimize: invalid bounds")
	// ErrBracket reports a scalar root bracket whose ends do not change sign.
	ErrBracket = errors.New("optimize: f(a) and f(b) must have different signs")
	// ErrInput reports any other invalid argument.
	ErrInput = errors.New("optimize: invalid input")
)

Errors returned for invalid input. Failure to converge is not an error: it is reported through the Success, Status and Message fields of a result, as SciPy does, so the caller still receives the best point found.

View Source
var ErrNotConverged = errors.New("optimize: optimal parameters not found")

ErrNotConverged is returned by CurveFit when the least-squares solver did not converge (SciPy raises RuntimeError("Optimal parameters not found")). The partial result is returned alongside it.

Functions

func ApproxFprime

func ApproxFprime(f func([]float64) float64, x []float64, epsilon float64) []float64

ApproxFprime returns the forward-difference gradient of f at x with the absolute step epsilon (0 selects SciPy's default sqrt(eps) ≈ 1.49e-8), as scipy.optimize.approx_fprime(xk, f, epsilon) does.

func Bracket

func Bracket(f func(float64) float64, xa, xb float64) (x [3]float64, fx [3]float64, nfev int, err error)

Bracket searches downhill from xa, xb (0 and 1 when both are zero) for three points xa, xb, xc with f(xb) below f(xa) and f(xc); it is scipy.optimize.bracket(func, xa, xb, grow_limit=110, maxiter=1000). It returns the points, their values and the number of evaluations, or a *BracketError.

func CheckGrad

func CheckGrad(f func([]float64) float64, grad func(x, g []float64), x []float64, epsilon float64) float64

CheckGrad returns the Euclidean norm of the difference between grad(x) and a forward-difference gradient of f with absolute step epsilon (0 for the default), like scipy.optimize.check_grad(func, grad, x0, epsilon=...). grad writes the gradient at x into its second argument.

Types

type Bound

type Bound struct {
	Lo, Hi float64
}

Bound is a closed interval [Lo, Hi] for one variable. Use math.Inf(-1) and math.Inf(1) for a missing side; the zero value Bound{} pins the variable to zero, so use Unbounded() for a free variable.

func Unbounded

func Unbounded() Bound

Unbounded returns the free interval (-Inf, +Inf).

type BracketError

type BracketError struct {
	Xa, Xb, Xc, Fa, Fb, Fc float64
	Nfev                   int
	Msg                    string
}

BracketError reports that Bracket found no valid bracket. It carries the last three points, as SciPy's BracketError does.

func (*BracketError) Error

func (e *BracketError) Error() string

type BracketOptions

type BracketOptions struct {
	Xtol, Rtol float64
	MaxIter    int
}

BracketOptions configures the bracketing root finders. Zero values select SciPy's defaults: Xtol 2e-12, Rtol 4·eps (the smallest allowed), MaxIter 100.

type BruteOptions

type BruteOptions struct {
	// Ns is the number of grid points per axis when Range.Step is 0 (20).
	Ns int
	// NoFinish skips the final local minimisation (SciPy's finish=None).
	NoFinish bool
	// Finish is the local method polishing the best grid point
	// (Nelder-Mead, SciPy's fmin). Finish options come from FinishOptions.
	Finish        Method
	FinishOptions *Options
	// Workers evaluates the grid with this many goroutines (1). The
	// objective must then be safe for concurrent use.
	Workers int
}

BruteOptions configures Brute.

type BruteResult

type BruteResult struct {
	// X and Fun are the minimiser and its value (after the finish step).
	X   []float64
	Fun float64
	// Success is false when the finish step did not converge.
	Success bool
	// Axes holds the coordinates of each grid axis; Jout the objective on
	// the grid, row-major over the axes (shape len(Axes[0]) × …).
	Axes [][]float64
	Jout []float64
}

BruteResult is Brute's full output (SciPy's full_output=True).

func Brute

func Brute(f func([]float64) float64, ranges []Range, o *BruteOptions) (*BruteResult, error)

Brute minimises f by evaluating it on a grid and polishing the best point with a local method; it is scipy.optimize.brute(func, ranges, Ns=20, full_output=True, finish=fmin, workers=...). The finish step may leave the grid's box, as in SciPy.

type CurveFitOptions

type CurveFitOptions struct {
	// Sigma is the uncertainty of ydata: a 1-D array of standard
	// deviations (the residuals are divided by it) or a 2-D covariance
	// matrix (the residuals are whitened by its Cholesky factor). nil means
	// unit weights.
	Sigma *ndarray.Array
	// AbsoluteSigma uses Sigma in an absolute sense: Pcov is not scaled by
	// the reduced chi-square.
	AbsoluteSigma bool
	// Bounds on the parameters; bounds select TRF.
	Bounds []Bound
	// Method is LM (the default without bounds) or TRF (with bounds).
	Method LSMethod
	// MaxFev, Ftol, Xtol and Gtol are passed to the solver (SciPy's maxfev
	// for LM, max_nfev for TRF; tolerances 1.49012e-8 for LM, 1e-8 for TRF).
	MaxFev           int
	Ftol, Xtol, Gtol float64
	// Loss and FScale select a robust loss (TRF only), as SciPy forwards
	// loss= and f_scale= to least_squares.
	Loss   Loss
	FScale float64
	// Workers evaluates the model over the samples with this many
	// goroutines when there are at least 2048 samples (default 1). The
	// model must then be safe for concurrent use; the result is identical.
	Workers int
}

CurveFitOptions configures CurveFit. Zero values select SciPy's defaults.

type CurveFitResult

type CurveFitResult struct {
	// Popt holds the fitted parameters.
	Popt []float64
	// Pcov is the n×n covariance of Popt, row-major. When it cannot be
	// estimated every entry is +Inf and CovOK is false (SciPy warns).
	Pcov  []float64
	CovOK bool
	// Fvec holds the (weighted) residuals at Popt and Nfev the solver's
	// count of residual evaluations (SciPy's infodict['nfev']: MINPACK's
	// count including finite differences for LM, least_squares' count
	// excluding them for TRF).
	Fvec []float64
	Nfev int
	// Status and Message are the solver's (SciPy's ier and mesg).
	Status  int
	Message string
}

CurveFitResult is the outcome of CurveFit: SciPy's (popt, pcov) plus its full_output information.

func CurveFit

func CurveFit(f func(x float64, p []float64) float64, xdata, ydata *ndarray.Array, p0 []float64, o *CurveFitOptions) (*CurveFitResult, error)

CurveFit fits a 1-D model f(x; p) to (xdata, ydata) by nonlinear least squares; it is scipy.optimize.curve_fit(f, xdata, ydata, p0, sigma=..., absolute_sigma=..., bounds=..., method=...). See CurveFitModel for models of several independent variables or with analytic derivatives.

Example

popt, pcov = scipy.optimize.curve_fit(lambda x, a, b, c: a*np.exp(-b*x) + c, x, y)

package main

import (
	"fmt"
	"math"

	"github.com/go-ndarray/ndarray"
	"github.com/go-numerics/optimize"
)

func main() {
	x, _ := ndarray.Linspace(0, 4, 50)
	y := x.Map(func(v float64) float64 { return 2.5*math.Exp(-1.3*v) + 0.5 + 0.01*math.Sin(37*v) })
	model := func(x float64, p []float64) float64 { return p[0]*math.Exp(-p[1]*x) + p[2] }
	r, err := optimize.CurveFit(model, x, y, []float64{1, 1, 1}, nil)
	if err != nil {
		panic(err)
	}
	perr := r.Perr()
	fmt.Printf("a = %.3f ± %.3f, b = %.3f ± %.3f, c = %.3f ± %.3f\n",
		r.Popt[0], perr[0], r.Popt[1], perr[1], r.Popt[2], perr[2])
}
Output:
a = 2.500 ± 0.004, b = 1.300 ± 0.005, c = 0.500 ± 0.002

func CurveFitModel

func CurveFitModel(model Model, xdata, ydata *ndarray.Array, p0 []float64, o *CurveFitOptions) (*CurveFitResult, error)

CurveFitModel is CurveFit for a general Model. xdata is 1-D (m samples) or 2-D (k, m); ydata is 1-D with m elements. p0 sets the number of parameters; nil starts from SciPy's default (ones, or a feasible point inside the bounds) when Bounds gives the count.

Without bounds the default method is MINPACK's Levenberg–Marquardt and the covariance is SciPy leastsq's cov_x; with TRF it is the pseudo-inverse of JᵀJ from the SVD of the final Jacobian, discarding singular values below eps·max(m, n)·s₀. Unless AbsoluteSigma is set the covariance is scaled by the reduced chi-square Σr²/(m − n).

func (*CurveFitResult) PcovArray

func (r *CurveFitResult) PcovArray() *ndarray.Array

PcovArray returns Pcov as an n×n ndarray.

func (*CurveFitResult) Perr

func (r *CurveFitResult) Perr() []float64

Perr returns the one-standard-deviation errors √diag(Pcov).

func (*CurveFitResult) PoptArray

func (r *CurveFitResult) PoptArray() *ndarray.Array

PoptArray returns Popt as a 1-D ndarray.

type FDMethod

type FDMethod int

FDMethod selects a finite-difference scheme.

const (
	// FD2Point is the forward difference (SciPy '2-point'), one evaluation
	// per variable, error O(h).
	FD2Point FDMethod = iota
	// FD3Point is the central difference (SciPy '3-point'), two evaluations
	// per variable, error O(h²); it falls back to a one-sided second-order
	// formula next to a bound.
	FD3Point
)

type LPConstraints

type LPConstraints struct {
	// AUb and BUb are the inequality constraints A_ub·x <= b_ub.
	AUb, BUb *ndarray.Array
	// AEq and BEq are the equality constraints A_eq·x == b_eq.
	AEq, BEq *ndarray.Array
	// Bounds gives one interval per variable. nil means linprog's default
	// (0, +Inf) for every variable — note that this differs from Minimize,
	// where nil means unbounded. Use Unbounded() for a free variable.
	Bounds []Bound
}

LPConstraints holds the constraints of a linear program

A_ub·x <= b_ub,  A_eq·x == b_eq,  Bounds[i].Lo <= x[i] <= Bounds[i].Hi,

in the shape of linprog's A_ub, b_ub, A_eq, b_eq and bounds arguments. Any of the matrices may be nil (no such constraints); a non-nil matrix must be 2-D with len(c) columns and its right-hand side 1-D with one entry per row.

type LPDual

type LPDual struct {
	// Residual is the slack of each constraint: b_ub − A_ub·x for
	// inequalities, b_eq − A_eq·x for equalities, x − lb and ub − x for
	// the bounds (±Inf for a missing bound).
	Residual []float64
	// Marginals is the sensitivity of the optimal objective to the
	// constraint's right-hand side (∂fun/∂b), with SciPy's HiGHS sign
	// convention: ineqlin and upper marginals are <= 0, lower marginals
	// are >= 0, eqlin marginals have any sign. They satisfy
	// c = A_ubᵀ·ineqlin + A_eqᵀ·eqlin + lower + upper. nil unless Status
	// is 0.
	Marginals []float64
}

LPDual pairs the residuals of one family of constraints with their marginals (dual values), like the ineqlin, eqlin, lower and upper fields of SciPy's HiGHS result.

type LPMethod

type LPMethod string

LPMethod names a linear-programming algorithm. The values are the strings SciPy's linprog accepts for the matching legacy methods.

const (
	// LPSimplex is a dense bounded-variable revised simplex method with an
	// LU-factorised basis and product-form updates. It is the default: it
	// ends on a vertex, its duals are exact for the final basis, and it
	// certifies infeasibility and unboundedness directly.
	LPSimplex LPMethod = "revised simplex"
	// LPInteriorPoint is the homogeneous self-dual primal–dual
	// interior-point method of Andersen & Andersen with Mehrotra's
	// predictor–corrector, solving dense normal equations by Cholesky. Its
	// solution is interior to the optimal face, not a vertex.
	LPInteriorPoint LPMethod = "interior-point"
)

Linear-programming methods implemented by Linprog.

type LPOptions

type LPOptions struct {
	// Method selects the algorithm; the default is LPSimplex.
	Method LPMethod
	// Tol is the primal and dual feasibility tolerance. Default 1e-9 for
	// LPSimplex and 1e-8 for LPInteriorPoint (SciPy's interior-point
	// default). As in SciPy, the final feasibility check uses 10·√Tol.
	Tol float64
	// MaxIter bounds the iterations (simplex pivots and bound flips over
	// both phases, or interior-point iterations). Default 5000 for
	// LPSimplex and 1000 for LPInteriorPoint, as in SciPy.
	MaxIter int
	// NoPresolve disables the presolve pass (linprog's presolve=False).
	NoPresolve bool
}

LPOptions configures Linprog. Zero values select the defaults listed with each field.

type LPResult

type LPResult struct {
	// X is the solution. When Status is not 0 it is the last point the
	// solver held (it may violate the constraints), or nil when presolve
	// stopped before any solve.
	X []float64
	// Fun is c·X.
	Fun float64
	// Slack is b_ub − A_ub·X and Con is b_eq − A_eq·X.
	Slack, Con []float64
	// IneqLin, EqLin, Lower and Upper hold residuals and marginals of the
	// inequality rows, equality rows, lower bounds and upper bounds.
	IneqLin, EqLin, Lower, Upper LPDual
	// Status is SciPy's code: 0 optimal, 1 iteration limit, 2 infeasible,
	// 3 unbounded, 4 numerical difficulties.
	Status int
	// Success reports Status == 0.
	Success bool
	// Message describes the reason for termination.
	Message string
	// Nit counts iterations over all phases.
	Nit int
}

LPResult is the outcome of Linprog, the counterpart of the OptimizeResult that scipy.optimize.linprog returns.

func Linprog

func Linprog(c []float64, cons LPConstraints, opts *LPOptions) (*LPResult, error)

Linprog minimises c·x subject to cons; it is the counterpart of scipy.optimize.linprog(c, A_ub, b_ub, A_eq, b_eq, bounds, method=...). SciPy's default HiGHS solvers are not available; opts.Method selects a dense revised simplex (default) or a dense interior-point method, after a light version of the presolve SciPy applies to its legacy methods: empty rows and columns and fixed variables are removed, trivial infeasibility and unboundedness are reported with SciPy's messages, and linearly dependent equality rows are dropped (or reported inconsistent). Unlike SciPy, singleton rows are left to the solver. The marginals are those of the original problem, with HiGHS's sign convention.

Invalid input (shapes, NaN or Inf coefficients, NaN bounds, unknown method) is an error wrapping ErrDimension, ErrInput, ErrBounds or ErrMethod. Infeasible, unbounded or unsolved problems are not errors: they are reported through Status, Success and Message.

Example

scipy.optimize.linprog(c, A_ub=A, b_ub=b, bounds=[(None, None), (-3, None)])

package main

import (
	"fmt"
	"math"

	"github.com/go-ndarray/ndarray"
	"github.com/go-numerics/optimize"
)

func main() {
	A, _ := ndarray.FromData([]float64{-3, 1, 1, 2}, 2, 2)
	b, _ := ndarray.FromData([]float64{6, 4}, 2)
	r, _ := optimize.Linprog([]float64{-1, 4}, optimize.LPConstraints{AUb: A, BUb: b,
		Bounds: []optimize.Bound{optimize.Unbounded(), {Lo: -3, Hi: math.Inf(1)}}}, nil)
	fmt.Printf("x = [%.1f %.1f], fun = %.1f, %s\n", r.X[0], r.X[1], r.Fun, r.Message)
}
Output:
x = [10.0 -3.0], fun = -22.0, Optimization terminated successfully.

func (*LPResult) XArray

func (r *LPResult) XArray() *ndarray.Array

XArray returns X as a 1-D ndarray (nil when X is nil).

type LSMethod

type LSMethod string

LSMethod names a least-squares algorithm; the values are SciPy's least_squares method strings.

const (
	// TRF is the trust-region reflective method (SciPy's default), which
	// supports bounds and robust loss functions.
	TRF LSMethod = "trf"
	// LM is the Levenberg–Marquardt method of MINPACK (lmdif/lmder); it
	// needs m >= n and supports neither bounds nor robust losses.
	LM LSMethod = "lm"
)

Least-squares methods.

type LSOptions

type LSOptions struct {
	// Bounds constrains the variables (TRF only). nil means unbounded.
	Bounds []Bound
	// Ftol, Xtol and Gtol are the termination tolerances (1e-8 each).
	// SciPy lets None disable one; here a negative value disables it.
	Ftol, Xtol, Gtol float64
	// MaxNfev bounds function evaluations (100·n; for LM it bounds
	// MINPACK's count, which leaves out finite-difference evaluations).
	MaxNfev int
	// XScale is the characteristic scale of each variable. nil means 1
	// (TRF) or 'jac' (LM, SciPy's default for that method); XScaleJac
	// selects SciPy's x_scale='jac'.
	XScale    []float64
	XScaleJac bool
	// Loss is the robust loss (TRF only; LossLinear by default) and FScale
	// its soft margin between inliers and outliers (1).
	Loss   Loss
	FScale float64
	// DiffStep is the relative finite-difference step (SciPy diff_step);
	// 0 selects the default of the scheme.
	DiffStep float64
	// FDMethod is the finite-difference scheme for TRF (FD2Point). LM uses
	// MINPACK's own forward differences.
	FDMethod FDMethod
}

LSOptions configures LeastSquares. Zero values select SciPy's defaults.

type LSProblem

type LSProblem struct {
	// Func writes the m residuals at x into f. Required.
	Func func(x, f []float64)
	// Jac writes the m×n Jacobian at x into jac, row-major. Optional:
	// finite differences are used when it is nil.
	Jac func(x, jac []float64)
	// M is the number of residuals. Required.
	M int
}

LSProblem is a nonlinear least-squares problem: minimise ½·Σ ρ(f_i(x)²) over x, with m residuals f_i and n variables.

type LSResult

type LSResult struct {
	// X is the solution.
	X []float64
	// Cost is ½·Σρ(f_i²) at X.
	Cost float64
	// Fun holds the m residuals at X.
	Fun []float64
	// Jac is the m×n Jacobian at X, row-major (modified by the loss, as in
	// SciPy).
	Jac []float64
	// Grad is the gradient of the cost at X.
	Grad []float64
	// Optimality is the first-order optimality measure (max-norm of the
	// scaled gradient, as SciPy defines it per method).
	Optimality float64
	// ActiveMask is -1 or +1 for a variable at its lower or upper bound,
	// 0 otherwise.
	ActiveMask []int
	// Nfev and Njev count residual and Jacobian evaluations.
	Nfev, Njev int
	// Status follows SciPy: -1 improper input, 0 max_nfev reached, 1 gtol,
	// 2 ftol, 3 xtol, 4 both ftol and xtol satisfied.
	Status int
	// Success is Status > 0.
	Success bool
	// Message describes the termination.
	Message string
}

LSResult is the outcome of LeastSquares, mirroring the OptimizeResult of scipy.optimize.least_squares.

func LeastSquares

func LeastSquares(p LSProblem, x0 []float64, method LSMethod, opts *LSOptions) (*LSResult, error)

LeastSquares solves a nonlinear least-squares problem from x0; it is scipy.optimize.least_squares(fun, x0, jac=..., bounds=..., method=..., ftol=..., xtol=..., gtol=..., x_scale=..., loss=..., f_scale=..., diff_step=..., max_nfev=...). method "" selects TRF.

Example

scipy.optimize.least_squares(fun, x0, bounds=([0, 0], [2, 0.5]), loss='soft_l1')

package main

import (
	"fmt"
	"math"

	"github.com/go-numerics/optimize"
)

func main() {
	t := []float64{0, 1, 2, 3, 4, 5}
	y := []float64{1.0, 0.62, 0.36, 0.23, 0.13, 3.0} // the last point is an outlier
	p := optimize.LSProblem{M: len(t), Func: func(x, f []float64) {
		for i := range t {
			f[i] = x[0]*math.Exp(-x[1]*t[i]) - y[i]
		}
	}}
	r, _ := optimize.LeastSquares(p, []float64{1, 0.1}, optimize.TRF, &optimize.LSOptions{
		Bounds: []optimize.Bound{{Lo: 0, Hi: 2}, {Lo: 0, Hi: 0.5}}, Loss: optimize.LossSoftL1, FScale: 0.1})
	fmt.Printf("x = [%.3f %.3f] active %v\n", r.X[0], r.X[1], r.ActiveMask)
}
Output:
x = [0.990 0.466] active [0 0]

func (*LSResult) JacArray

func (r *LSResult) JacArray() *ndarray.Array

JacArray returns Jac as an m×n ndarray.

func (*LSResult) XArray

func (r *LSResult) XArray() *ndarray.Array

XArray returns X as a 1-D ndarray.

type LeastsqOptions

type LeastsqOptions struct {
	// Ftol is the relative error desired in the sum of squares and Xtol the
	// relative error desired in the solution (1.49012e-8 each; a negative
	// value means exactly 0). Gtol is the orthogonality desired between the
	// residuals and the Jacobian columns (0).
	Ftol, Xtol, Gtol float64
	// Maxfev bounds the calls to Func: 0 selects 100·(n+1) with a Jacobian
	// and 200·(n+1) without.
	Maxfev int
	// Epsfcn sets the forward-difference step sqrt(Epsfcn)·|x_j| when there
	// is no Jacobian; 0 selects the machine epsilon.
	Epsfcn float64
	// Factor sets the initial step bound Factor·|diag·x0| (100).
	Factor float64
	// Diag holds n positive scale factors for the variables; nil lets
	// MINPACK scale by the Jacobian column norms.
	Diag []float64
}

LeastsqOptions configures Leastsq. Zero values select the defaults of scipy.optimize.leastsq.

type LeastsqResult

type LeastsqResult struct {
	// X is the solution (or the last iterate).
	X []float64
	// CovX is inv(JᵀJ) at X, n×n row-major, from the final R and
	// permutation; nil when R is singular or the run did not converge
	// (Ier not in 1..4). Multiply it by the residual variance to obtain the
	// covariance of X.
	CovX []float64
	// Nfev counts calls to Func, including those of the finite-difference
	// Jacobian; Njev counts calls to Jac (0 without one).
	Nfev, Njev int
	// Fvec holds the residuals at X.
	Fvec []float64
	// Fjac is SciPy's infodict['fjac']: an n×m row-major array whose rows
	// are the columns of MINPACK's final QR workspace, so that the upper
	// triangle of its transpose's first n rows is R with Jac·P = Q·R.
	Fjac []float64
	// Ipvt is the permutation P: column j of P is column Ipvt[j]
	// (zero-based) of the identity.
	Ipvt []int
	// Qtf is the first n components of Qᵀ·Fvec.
	Qtf []float64
	// Ier is MINPACK's info: 1 to 4 mean a solution was found, 5 to 8 that
	// it was not.
	Ier int
	// Mesg is SciPy's message for Ier.
	Mesg string
}

LeastsqResult is the full output of Leastsq: SciPy's tuple (x, cov_x, infodict, mesg, ier).

func Leastsq

func Leastsq(p LSProblem, x0 []float64, opts *LeastsqOptions) (*LeastsqResult, error)

Leastsq minimises the sum of squares of p.M residuals with MINPACK's Levenberg–Marquardt method; it is scipy.optimize.leastsq(func, x0, Dfun=p.Jac, full_output=True, ftol, xtol, gtol, maxfev, epsfcn, factor, diag). Without p.Jac it runs lmdif (forward differences), otherwise lmder. Non-convergence is reported through Ier and Mesg.

func (*LeastsqResult) CovXArray

func (r *LeastsqResult) CovXArray() *ndarray.Array

CovXArray returns CovX as an n×n ndarray, or nil when absent.

func (*LeastsqResult) FjacArray

func (r *LeastsqResult) FjacArray() *ndarray.Array

FjacArray returns Fjac as an n×m ndarray.

func (*LeastsqResult) Success

func (r *LeastsqResult) Success() bool

Success reports whether Ier is 1, 2, 3 or 4, as SciPy's leastsq judges it.

func (*LeastsqResult) XArray

func (r *LeastsqResult) XArray() *ndarray.Array

XArray returns X as a 1-D ndarray.

type Loss

type Loss string

Loss names a robust loss ρ for LeastSquares (SciPy's loss argument).

const (
	LossLinear Loss = "linear"  // ρ(z) = z
	LossSoftL1 Loss = "soft_l1" // ρ(z) = 2((1+z)^½ − 1)
	LossHuber  Loss = "huber"   // ρ(z) = z if z <= 1 else 2z^½ − 1
	LossCauchy Loss = "cauchy"  // ρ(z) = ln(1+z)
	LossArctan Loss = "arctan"  // ρ(z) = arctan(z)
)

Loss functions, with z = (f/FScale)².

type Method

type Method string

Method names a minimisation algorithm. The values are the strings SciPy's minimize accepts for its method argument.

const (
	NelderMead Method = "Nelder-Mead"
	Powell     Method = "Powell"
	CG         Method = "CG"
	BFGS       Method = "BFGS"
	NewtonCG   Method = "Newton-CG"
	LBFGSB     Method = "L-BFGS-B"
)

Methods implemented by Minimize.

type Model

type Model struct {
	// F returns the model value at sample x for parameters p. Required.
	F func(x, p []float64) float64
	// Grad writes ∂f/∂p at sample x into grad (len(p) values). Optional:
	// finite differences are used when it is nil (SciPy's jac=callable).
	Grad func(x, p, grad []float64)
	// contains filtered or unexported fields
}

Model is a curve-fitting model y = f(x; p) evaluated at one sample. For 1-D xdata x has one element; for xdata of shape (k, m) — k independent variables, m samples, SciPy's layout — x holds the k values of one sample.

type NewtonOptions

type NewtonOptions struct {
	// Fprime is the derivative; nil selects the secant method.
	Fprime func(float64) float64
	// Fprime2 is the second derivative; with Fprime it selects Halley's
	// method.
	Fprime2 func(float64) float64
	// X1 is the second starting point of the secant method; 0 means
	// x0·(1+1e-4) ± 1e-4, as in SciPy. Set X1Set to use an explicit 0.
	X1    float64
	X1Set bool
	// Tol and Rtol are the absolute and relative tolerances on successive
	// iterates (numpy.isclose semantics).
	Tol, Rtol float64
	MaxIter   int
}

NewtonOptions configures Newton. Zero values select SciPy's defaults: Tol 1.48e-8, Rtol 0, MaxIter 50.

type Options

type Options struct {
	// Bounds constrains the variables (L-BFGS-B, Nelder-Mead and Powell, as in
	// SciPy). nil means unbounded.
	Bounds []Bound
	// Tol is minimize's tol argument: it sets the method's main tolerance
	// the way SciPy does (Nelder-Mead: Xatol and Fatol; Powell: Xtol and
	// Ftol; CG and BFGS: Gtol; Newton-CG: Xtol; L-BFGS-B: Ftol and Gtol)
	// when those fields are zero.
	Tol float64
	// MaxIter bounds the iterations (Nelder-Mead 200·n, Powell 1000·n,
	// CG, BFGS and Newton-CG 200·n, L-BFGS-B 15000). For Nelder-Mead and
	// Powell, setting only one of MaxIter and MaxFev leaves the other
	// unlimited, as in SciPy.
	MaxIter int
	// MaxFev bounds function evaluations (Nelder-Mead 200·n, Powell 1000·n,
	// L-BFGS-B 15000).
	MaxFev int
	// Gtol is the gradient tolerance: CG and BFGS stop when the Norm-norm
	// of the gradient is below it (1e-5); L-BFGS-B stops when the max-norm
	// of the projected gradient is below it (SciPy pgtol, 1e-5).
	Gtol float64
	// Ftol is the relative reduction tolerance: Powell (1e-4),
	// L-BFGS-B ((f_k - f_{k+1})/max(|f_k|,|f_{k+1}|,1) <= Ftol, 2.22e-9).
	Ftol float64
	// Xtol is the step tolerance: Powell (1e-4), Newton-CG (1e-5).
	Xtol float64
	// Xatol and Fatol are Nelder-Mead's absolute tolerances (1e-4 each).
	Xatol, Fatol float64
	// Adaptive makes Nelder-Mead adapt its coefficients to the dimension.
	Adaptive bool
	// InitialSimplex is Nelder-Mead's initial simplex, n+1 points.
	InitialSimplex [][]float64
	// Direc is Powell's initial set of directions, n vectors.
	Direc [][]float64
	// Norm is the order of the gradient norm for CG and BFGS (+Inf).
	Norm float64
	// Eps is the absolute finite-difference step for gradient-based
	// methods when no gradient is given (sqrt(eps) ≈ 1.49e-8 for CG and
	// BFGS, 1e-8 for L-BFGS-B). It is ignored when FiniteDiffRelStep is set.
	Eps float64
	// FiniteDiffRelStep selects a relative finite-difference step instead
	// of the absolute Eps (SciPy finite_diff_rel_step with jac='2-point').
	FiniteDiffRelStep float64
	// FDMethod is the finite-difference scheme (FD2Point).
	FDMethod FDMethod
	// M is the number of L-BFGS-B correction pairs (10).
	M int
	// MaxLS bounds L-BFGS-B line-search steps per iteration (20).
	MaxLS int
	// C1 and C2 are the Wolfe constants for CG (1e-4, 0.4) and BFGS
	// (1e-4, 0.9).
	C1, C2 float64
	// HessInv0 is BFGS's initial inverse Hessian, n*n row-major (identity).
	HessInv0 []float64
	// XRtol is BFGS's relative step tolerance (0: disabled).
	XRtol float64
	// Callback is called after each iteration with a copy of the current
	// point and the objective value; returning true stops the minimisation,
	// which then reports Success false and Status 99 (SciPy's callback
	// raising StopIteration).
	Callback func(x []float64, f float64) (stop bool)
}

Options configures Minimize. Zero values select the SciPy defaults of the chosen method, listed with each field.

type Problem

type Problem struct {
	// Func returns f(x).
	Func func(x []float64) float64
	// Grad writes ∇f(x) into grad (SciPy's jac=callable).
	Grad func(x, grad []float64)
	// FuncGrad returns f(x) and writes ∇f(x) into grad in one call
	// (SciPy's jac=True). It takes precedence over Func and Grad.
	FuncGrad func(x, grad []float64) float64
	// Hess writes the n×n Hessian at x into hess, row-major (SciPy's hess).
	Hess func(x, hess []float64)
	// HessP writes the Hessian at x times p into hp (SciPy's hessp).
	HessP func(x, p, hp []float64)
}

Problem describes a scalar objective for Minimize. Func is required unless FuncGrad is set. Derivatives are optional: when neither Grad nor FuncGrad is given, gradient-based methods use finite differences with SciPy's step rules (see Options.Eps, Options.FiniteDiffRelStep, Options.FDMethod).

Every callback receives the point first and writes its outputs into the slices that follow; it must not retain x.

type Range

type Range struct {
	Lo, Hi, Step float64
}

Range is one axis of a Brute grid. With Step 0 the axis holds Ns points from Lo to Hi inclusive (SciPy's (low, high) tuple, i.e. slice(low, high, Ns*1j)); with Step > 0 it holds Lo, Lo+Step, … below Hi (slice(low, high, step)).

type Result

type Result struct {
	// X is the solution.
	X []float64
	// Fun is the objective value at X.
	Fun float64
	// Jac is the gradient at X, when the method computes one.
	Jac []float64
	// HessInv is the final inverse-Hessian approximation as an n*n row-major
	// matrix (BFGS only, nil otherwise).
	HessInv []float64
	// Nit, Nfev, Njev and Nhev count iterations and evaluations of the
	// objective, the gradient and the Hessian. Function evaluations spent on
	// finite-difference gradients are counted in Nfev, and each such gradient
	// also counts once in Njev, as in SciPy.
	Nit, Nfev, Njev, Nhev int
	// Status is a method-specific termination code; 0 means success.
	Status int
	// Success reports whether the method met its convergence criteria.
	Success bool
	// Message describes the reason for termination.
	Message string
}

Result is the outcome of Minimize, the counterpart of SciPy's OptimizeResult for scipy.optimize.minimize.

func Minimize

func Minimize(p Problem, x0 []float64, method Method, opts *Options) (*Result, error)

Minimize minimises p from x0 with the given method; it is scipy.optimize.minimize(fun, x0, method=..., jac=..., hess=..., bounds=..., tol=..., options=...).

Invalid input returns an error; non-convergence returns a Result with Success false.

Example

scipy.optimize.minimize(rosen, x0, jac=rosen_der, method='BFGS')

package main

import (
	"fmt"
	"math"

	"github.com/go-numerics/optimize"
)

func main() {
	rosen := func(x []float64) float64 {
		return 100*math.Pow(x[1]-x[0]*x[0], 2) + math.Pow(1-x[0], 2)
	}
	grad := func(x, g []float64) {
		g[0] = -400*x[0]*(x[1]-x[0]*x[0]) - 2*(1-x[0])
		g[1] = 200 * (x[1] - x[0]*x[0])
	}
	r, err := optimize.Minimize(optimize.Problem{Func: rosen, Grad: grad}, []float64{-1.2, 1}, optimize.BFGS, nil)
	if err != nil {
		panic(err)
	}
	fmt.Printf("x = [%.6f %.6f], nit = %d, %s\n", r.X[0], r.X[1], r.Nit, r.Message)
}
Output:
x = [1.000000 1.000000], nit = 32, Optimization terminated successfully.
Example (Bounds)

scipy.optimize.minimize(f, x0, method='L-BFGS-B', bounds=[(-2, 0.5), (None, None)])

package main

import (
	"fmt"

	"github.com/go-numerics/optimize"
)

func main() {
	f := func(x []float64) float64 { return (x[0]-1)*(x[0]-1) + (x[1]+2)*(x[1]+2) }
	opts := &optimize.Options{Bounds: []optimize.Bound{{Lo: -2, Hi: 0.5}, optimize.Unbounded()}}
	r, _ := optimize.Minimize(optimize.Problem{Func: f}, []float64{0, 0}, optimize.LBFGSB, opts)
	fmt.Printf("x = [%.4f %.4f]\n", r.X[0], r.X[1])
}
Output:
x = [0.5000 -2.0000]

func (*Result) HessInvArray

func (r *Result) HessInvArray() *ndarray.Array

HessInvArray returns HessInv as an n×n ndarray, or nil when absent.

func (*Result) XArray

func (r *Result) XArray() *ndarray.Array

XArray returns X as a 1-D ndarray.

type RootOptions

type RootOptions struct {
	// Xtol stops the iteration when the relative error between two
	// consecutive iterates is at most Xtol (1.49012e-8; a negative value
	// means exactly 0).
	Xtol float64
	// Maxfev bounds the calls to Func: 0 selects 200·(n+1) without a
	// Jacobian (forward-difference calls included) and 100·(n+1) with one.
	Maxfev int
	// Eps sets the forward-difference step sqrt(Eps)·|x_j| when there is no
	// Jacobian (SciPy's eps for root, epsfcn for fsolve); 0 selects the
	// machine epsilon.
	Eps float64
	// Factor sets the initial step bound Factor·|Diag·x0| (100).
	Factor float64
	// Diag holds n positive scale factors for the variables; nil lets
	// MINPACK scale by the Jacobian column norms.
	Diag []float64
}

RootOptions configures Root and Fsolve. Zero values select SciPy's defaults.

type RootResult

type RootResult struct {
	// X is the solution (or the last iterate).
	X []float64
	// Fun holds F(X).
	Fun []float64
	// Fjac is SciPy's fjac: the orthogonal factor Q of the QR factorisation
	// of the final approximate Jacobian, stored column-wise, so that as an
	// n×n row-major array it reads Qᵀ (row i is column i of Q).
	Fjac []float64
	// R is the upper-triangular factor packed by rows (n(n+1)/2 entries).
	R []float64
	// Qtf is Qᵀ·Fun.
	Qtf []float64
	// Nfev counts calls to Func, including those of finite differences;
	// Njev counts calls to Jac (0 without one). SciPy's root and fsolve
	// report two more calls: those their wrappers make outside MINPACK (a
	// shape check and the initial evaluation), which do not happen here.
	Nfev, Njev int
	// Status is MINPACK's info: 1 converged, 2 Maxfev reached, 3 Xtol too
	// small, 4 and 5 no progress.
	Status int
	// Success is Status == 1.
	Success bool
	// Message is SciPy's message for Status.
	Message string
}

RootResult is the outcome of Root and Fsolve, mirroring the OptimizeResult of scipy.optimize.root(method='hybr').

func Fsolve

func Fsolve(p SystemProblem, x0 []float64, opts *RootOptions) (*RootResult, error)

Fsolve is scipy.optimize.fsolve(func, x0, fprime=p.Jac, full_output=True, xtol, maxfev, epsfcn=opts.Eps, factor, diag): the same computation as Root, returning the same result (SciPy's infodict holds Fun as fvec, Nfev, Njev, Fjac, R and Qtf; ier is Status and mesg Message).

func Root

func Root(p SystemProblem, x0 []float64, opts *RootOptions) (*RootResult, error)

Root solves the square system p.Func(x) = 0 from x0 with Powell's hybrid method of MINPACK; it is scipy.optimize.root(fun, x0, jac=p.Jac, method='hybr', options={'xtol', 'maxfev', 'eps', 'factor', 'diag'}) (the 'band' and 'col_deriv' options are not supported). Without p.Jac it runs hybrd, otherwise hybrj. Non-convergence is reported through Status and Message.

Example

scipy.optimize.root(fun, [0, 0], method='hybr')

package main

import (
	"fmt"
	"math"

	"github.com/go-numerics/optimize"
)

func main() {
	fun := func(x, f []float64) {
		f[0] = x[0] + 0.5*math.Pow(x[0]-x[1], 3) - 1
		f[1] = 0.5*math.Pow(x[1]-x[0], 3) + x[1]
	}
	r, _ := optimize.Root(optimize.SystemProblem{Func: fun}, []float64{0, 0}, nil)
	fmt.Printf("x = [%.6f %.6f] %v\n", r.X[0], r.X[1], r.Success)
}
Output:
x = [0.841164 0.158836] true

func (*RootResult) FjacArray

func (r *RootResult) FjacArray() *ndarray.Array

FjacArray returns Fjac as an n×n ndarray.

func (*RootResult) XArray

func (r *RootResult) XArray() *ndarray.Array

XArray returns X as a 1-D ndarray.

type RootScalarOptions

type RootScalarOptions struct {
	// Bracket [a, b] for the bracketing methods.
	Bracket []float64
	// Fprime and Fprime2 are the first and second derivatives.
	Fprime, Fprime2 func(float64) float64
	// X0 and X1 are starting points for Newton, Halley and the secant
	// method; X0Set/X1Set mark them as given (0 is a valid start).
	X0, X1       float64
	X0Set, X1Set bool
	// Xtol, Rtol and MaxIter override the method's defaults.
	Xtol, Rtol float64
	MaxIter    int
}

RootScalarOptions configures RootScalar.

type RootScalarResult

type RootScalarResult struct {
	// Root is the estimated root.
	Root float64
	// Iterations and FunctionCalls count the work done (derivative
	// evaluations count as function calls for Newton and Halley).
	Iterations, FunctionCalls int
	// Converged reports whether the tolerance was met.
	Converged bool
	// Flag is "converged", "convergence error" or a message describing the
	// failure.
	Flag string
	// Method is the algorithm used.
	Method string
}

RootScalarResult describes a scalar root search, like SciPy's RootResults.

func Bisect

func Bisect(f func(float64) float64, a, b float64, o *BracketOptions) (*RootScalarResult, error)

Bisect finds a root of f in [a, b] by bisection; it is scipy.optimize.bisect.

func Brenth

func Brenth(f func(float64) float64, a, b float64, o *BracketOptions) (*RootScalarResult, error)

Brenth is Brentq with hyperbolic instead of inverse quadratic extrapolation (Bus & Dekker 1975); it is scipy.optimize.brenth.

func Brentq

func Brentq(f func(float64) float64, a, b float64, o *BracketOptions) (*RootScalarResult, error)

Brentq finds a root of f in [a, b], where f(a) and f(b) have opposite signs, with Brent's (1973) method using inverse quadratic extrapolation; it is scipy.optimize.brentq(f, a, b, xtol, rtol, maxiter, full_output=True). The returned root x satisfies |x − x₀| <= 4·(Xtol + Rtol·|x₀|) for the exact root x₀.

A bracket without a sign change returns ErrBracket; a NaN function value returns an error together with the result describing where it happened; running out of iterations is reported by Converged false.

Example

scipy.optimize.brentq(lambda x: x**3 - 2*x - 5, 2, 3)

package main

import (
	"fmt"

	"github.com/go-numerics/optimize"
)

func main() {
	r, _ := optimize.Brentq(func(x float64) float64 { return x*x*x - 2*x - 5 }, 2, 3, nil)
	fmt.Printf("%.12f after %d iterations\n", r.Root, r.Iterations)
}
Output:
2.094551481542 after 7 iterations

func Newton

func Newton(f func(float64) float64, x0 float64, o *NewtonOptions) (*RootScalarResult, error)

Newton finds a root of f near x0 with the Newton–Raphson method (given Fprime), Halley's method (given Fprime and Fprime2) or the secant method (neither); it is scipy.optimize.newton(func, x0, fprime, tol=1.48e-8, maxiter=50, fprime2, x1, rtol=0, full_output=True, disp=False).

Failure to converge (including a zero derivative) is reported by Converged false and Flag "convergence error", not as an error.

func Ridder

func Ridder(f func(float64) float64, a, b float64, o *BracketOptions) (*RootScalarResult, error)

Ridder finds a root of f in [a, b] with Ridders' (1979) method; it is scipy.optimize.ridder.

func RootScalar

func RootScalar(f func(float64) float64, method string, o *RootScalarOptions) (*RootScalarResult, error)

RootScalar finds a root of a scalar function; it is scipy.optimize.root_scalar(f, method, bracket, fprime, fprime2, x0, x1, xtol, rtol, maxiter). method is one of "brentq", "brenth", "ridder", "bisect", "newton", "halley", "secant", or "" to choose as SciPy does: a bracket selects brentq; otherwise x0 with Fprime selects newton (halley with Fprime2), x0 with x1 the secant method, and x0 alone Newton with a finite-difference derivative.

type ScalarMethod

type ScalarMethod string

ScalarMethod names a scalar minimisation algorithm; the values are the method strings of scipy.optimize.minimize_scalar.

const (
	// Brent is Brent's method: parabolic interpolation safeguarded by
	// golden-section steps, from a bracket.
	Brent ScalarMethod = "brent"
	// Bounded is Brent's method restricted to a finite interval (the
	// algorithm of fminbound).
	Bounded ScalarMethod = "bounded"
	// Golden is the golden-section search, from a bracket.
	Golden ScalarMethod = "golden"
)

Scalar minimisation methods.

type ScalarOptions

type ScalarOptions struct {
	// Bracket is a bracketing interval for Brent and Golden: two points
	// (xa, xb) from which Bracket searches downhill, or three points
	// (xa, xb, xc) with xa < xb < xc and f(xb) < f(xa), f(xc). nil starts
	// the downhill search from (0, 1).
	Bracket []float64
	// Bounds is the interval for the Bounded method; setting it selects
	// Bounded when the method is "".
	Bounds *Bound
	// Tol is minimize_scalar's tol: it sets Xtol (Brent, Golden) or Xatol
	// (Bounded) when those are zero.
	Tol float64
	// Xtol is the relative tolerance of Brent (1.48e-8) and Golden
	// (√eps ≈ 1.49e-8).
	Xtol float64
	// Xatol is the absolute tolerance of Bounded (1e-5).
	Xatol float64
	// MaxIter bounds iterations (Brent 500, Golden 5000) or function
	// evaluations (Bounded 500).
	MaxIter int
}

ScalarOptions configures MinimizeScalar. Zero values select SciPy's defaults.

type ScalarResult

type ScalarResult struct {
	X, Fun    float64
	Nit, Nfev int
	Status    int // Bounded only: 0 solution found, 1 maxfun, 2 NaN
	Success   bool
	Message   string
}

ScalarResult is the outcome of MinimizeScalar, the OptimizeResult of scipy.optimize.minimize_scalar.

func Fminbound

func Fminbound(f func(float64) float64, x1, x2, xtol float64, maxfun int) (*ScalarResult, error)

Fminbound minimises f on [x1, x2] with absolute tolerance xtol (1e-5 when 0) and at most maxfun evaluations (500 when 0); it is scipy.optimize.fminbound(func, x1, x2, xtol, maxfun, full_output=True).

func MinimizeScalar

func MinimizeScalar(f func(float64) float64, method ScalarMethod, o *ScalarOptions) (*ScalarResult, error)

MinimizeScalar minimises a function of one variable; it is scipy.optimize.minimize_scalar(fun, bracket=..., bounds=..., method=..., tol=..., options={'xtol'/'xatol', 'maxiter'}). method "" selects Brent, or Bounded when o.Bounds is set.

A failed bracket search is not an error: as in SciPy, the result then holds the best of the three points and Success is false.

type SystemProblem

type SystemProblem struct {
	// Func writes the n values of F at x into f. Required. It must not
	// modify x.
	Func func(x, f []float64)
	// Jac writes the n×n Jacobian of F at x into jac, row-major
	// (jac[i*n+j] = ∂F_i/∂x_j). Optional: MINPACK's forward differences are
	// used when it is nil.
	Jac func(x, jac []float64)
}

SystemProblem is a square system of nonlinear equations F(x) = 0 with n equations in n unknowns.

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL