SF_OPTIMIZE

Description

SciFortran module for minimization and root finding

Quick access

Routines:

adaptive_mix(), bisect(), bracket(), brent(), brentq(), broyden1(), broyden_mix(), curvefit(), dbrent(), fmin(), fmin_bfgs(), fmin_cg(), fmin_cgminimize(), fmin_cgplus(), fsolve(), fzero(), leastsq(), linear_mix(), newton()

Used modules

  • cgfit_routines

  • broyden_routines

  • sf_constants: contains mathematical, numerical and physical constants

  • sf_linalg: SciFortran module for linear algebra

Subroutines and functions

interface  sf_optimize/fmin_cg(p, f, df, iter, fret[, ftol, itmax, istop, iverbose, err, deps])

Generic interface to minimize a function of several variables with the Fletcher-Reeves-Polak-Ribiere conjugate gradient method (Numerical Recipes 10.6), with the line minimizations done with the Brent method. The gradient is supplied by the user (fmin_cg_df) or computed numerically (fmin_cg_f). The location of the minimum is returned in p, the number of iterations in iter and the value of the function in fret.

The convergence is controlled by istop and ftol: with \(a = |F_n-F_{n-1}|/(1+|F_n|)\) and \(b = \|p_n-p_{n-1}\|^2/(1+\|p_n\|^2)\), the iteration stops when \(a<ftol\) and \(b<ftol\) (istop=0), when \(a<ftol\) (istop=1) or when \(b<ftol\) (istop=2). The last values of \(a\) and \(b\) are returned in err.

Parameters:

  • p (•) [real, inout] – in: starting point; out: location of the minimum

  • f [real] – function to be minimized, real(8) function f(x) with real(8) x(:)

  • df [real] – gradient of f, real(8) function df(x) returning a vector with the size of x

  • iter [integer, out] – out: number of iterations performed

  • fret [real, out] – out: value of f at the minimum

Options:

  • ftol [real] – convergence tolerance, default 1d-5

  • itmax [integer] – maximum number of iterations, default 500

  • istop [integer] – stopping criterion, 0 (default) both a and b < ftol, 1 a < ftol, 2 b < ftol

  • iverbose [logical] – print information on the convergence, default .false.

  • err (2) [real] – out: last values of the convergence measures [a,b]

  • deps [real] – step of the numerical derivative, default 0

interface  sf_optimize/fmin_cgplus(p, func, fjac, iter, fret, fcn[, ftol, itmax, imethod, iverb1, iverb2, deps])

Generic interface to minimize a function of several variables with the nonlinear conjugate gradient method of the CG+ package (Gilbert and Nocedal), a very old Fortran 77 code. The gradient is supplied by the user (fmin_cgplus_df) or computed numerically (fmin_cgplus_f). The variant of the method is selected by imethod, and the iteration stops when every component of the gradient is smaller than \(ftol(1+|f|)\). The location of the minimum is returned in p.

Parameters:

  • p (•) [real, inout] – in: starting point; out: location of the minimum

  • func [external] – function to be minimized, real(8) function func(x) with real(8) x(:)

  • fjac [external] – gradient of func, real(8) function fjac(x) returning a vector with the size of x

  • iter [integer, out] – out: number of iterations performed

  • fret [real] – out: value of func at the minimum

  • fcn [real] – function to be minimized, real(8) function fcn(x) with real(8) x(:)

Options:

  • ftol [real] – convergence tolerance on the gradient, default 1d-5

  • itmax [integer] – maximum number of iterations, default 1000

  • imethod [integer] – 1 Fletcher-Reeves, 2 Polak-Ribiere (default), 3 positive Polak-Ribiere

  • iverb1 [integer] – output frequency, < 0 none (default), 0 first and last iteration, > 0 every iverb1 iterations

  • iverb2 [integer] – amount of output, 0 (default) to 3: 1 initial x and gradient, 2 x every iteration, 3 x and gradient

  • deps [real] – optional: step of the numerical derivative, default 0

interface  sf_optimize/fmin_cgminimize(p, fcn, iter, fret[, ftol, itmax, iverbose, mode, new_version, hh_par])

Generic interface to minimize a function of several variables with a conjugate gradient routine adapted from a very old Fortran 77 code, a quasi-Newton minimizer with numerical gradient. The function is a function (fmin_cgminimize_func) or a subroutine (fmin_cgminimize_sub). Two versions of the minimizer are available, selected by new_version. The location of the minimum is returned in p.

Parameters:

  • p (•) [real, inout] – in: starting point; out: location of the minimum

  • fcn [real] – function to be minimized, real(8) function fcn(x) with real(8) x(:)

  • iter [integer] – out: number of iterations performed

  • fret [real] – out: value of the function at the minimum

Options:

  • ftol [real] – convergence tolerance, default 1d-5

  • itmax [integer] – maximum number of iterations, default 1000 (calls: 2*N*itmax)

  • iverbose [logical] – print information, default .false.

  • mode [integer] – mode of the minimizer, default 1 (the value passed is currently ignored)

  • new_version [logical] – .true. uses the new minimizer, .false. (default) the old one

  • hh_par [real] – step of the numerical gradient, default 1d-5

interface  sf_optimize/leastsq(func, a, m, dfunc[, tol, info])

Generic interface to minimize the sum of squares of \(m\) nonlinear functions of \(n\) unknowns, \(m\ge n\), i.e. to find a so that \(f(a)^T f(a)\) is minimum, with the Levenberg-Marquardt routines of MINPACK. The Jacobian is computed numerically (lmdif, leastsq_lmdif_func and leastsq_lmdif_sub) or supplied by the user (lmder, leastsq_lmder_func and leastsq_lmder_sub). The functions are a function func(a,m) returning a vector of size \(m\), or a subroutine func(a,m,f); the Jacobian dfunc(a,m) is a matrix of shape \([m,n]\).

The program stops, writing the reason to the file LEASTSQ_ERROR.err, if the status returned by MINPACK is not 1.

Parameters:

  • func [external]

  • a (•) [real] – in: starting guess; out: solution, size n

  • m [integer] – number of functions, m >= n

  • dfunc [external]

Options:

  • tol [real] – optional: tolerance, default 1d-15

  • info [integer] – optional out: status of the MINPACK routine

interface  sf_optimize/curvefit(model_func, a, xdata, ydata, model_dfunc[, tol, info])

Generic interface to fit a model function to data with the non-linear least squares Levenberg-Marquardt routines of MINPACK: the parameters a of model_func(x,a) are adjusted to minimize the sum of the squares of model_func(xdata,a)-ydata. The Jacobian is computed numerically (lmdif, curvefit_lmdif_func and curvefit_lmdif_sub) or supplied by the user (lmder, curvefit_lmder_func and curvefit_lmder_sub). The model is a function model_func(x,a) returning a vector with the size of x, or a subroutine model_func(x,a,f); the derivative model_dfunc(x,a) is a matrix of shape \([size(x),size(a)]\). The model is evaluated on the whole array xdata at once.

Parameters:

  • model_func [external]

  • a (•) [real] – in: starting guess of the parameters; out: fitted parameters, size n

  • xdata (•) [real, in,required] – abscissas of the data, size m

  • ydata (size(xdata)) [real] – ordinates of the data, size m

  • model_dfunc [external]

Options:

  • tol [real] – optional: tolerance, default 1d-15

  • info [integer] – optional out: status of the MINPACK routine

interface  sf_optimize/dbrent(func, dfunc, xmin[, brack, tol, niter])

Generic interface to find the minimum of a function of one variable, func, with the Brent method using the derivative: the user supplies the derivative dfunc (dbrent_wgrad), or it is computed numerically (dbrent_nograd). The abscissa of the minimum is returned in xmin. The bracketing of the minimum, brack, works as in brent: a bracketing triplet, or two starting points from which it is searched with bracket, with default points 0 and 1. The program stops if brack has a single element.

Parameters:

  • func [external]

  • dfunc [external]

  • xmin [real, inout] – abscissa of the minimum

Options:

  • brack (•) [real] – optional: 2 starting points or 3 bracketing points, default [0,1]

  • tol [real] – optional: fractional precision of the minimum, default 1d-9

  • niter [integer] – optional: maximum number of iterations, default 200

interface  sf_optimize/fmin_bfgs(func, grad, x[, l, u, nbd, factr, pgtol, iprint, nloop])

Generic interface to minimize a function of several variables with the limited-memory BFGS algorithm with bounds, L-BFGS-B (Zhu, Byrd, Lu and Nocedal), subject to the optional bounds \(l\le x\le u\). The gradient is supplied by the user (bfgs_with_grad) or computed numerically (bfgs_no_grad). The type of the bounds is given by nbd, the bounds are ignored if nbd is absent. The location of the minimum is returned in x.

Parameters:

  • func [external]

  • grad [external]

  • x (•) [real, allocatable] – in: starting point; out: location of the minimum

Options:

  • l (•) [real, allocatable] – optional: lower bounds, size n

  • u (•) [real, allocatable] – optional: upper bounds, size n

  • nbd (•) [integer, allocatable] – optional: 0 none, 1 lower, 2 both, 3 upper bound

  • factr [real] – optional: f tolerance in eps units, default 1d7

  • pgtol [real] – optional: projected gradient tolerance, default 1d-5

  • iprint [integer] – optional: <0 (default) no output, larger prints more

  • nloop [integer] – optional: maximum number of loops

interface  sf_optimize/linear_mix(x, fx, alpha)

Generic interface for the linear mixing of an array of rank 1 to 7, real (d_linear_mix_1 to d_linear_mix_7) or complex (c_linear_mix_1 to c_linear_mix_7): \(x \leftarrow x + \alpha F_x\), with \(F_x = x_{out}-x_{in}\) the residual of the iteration.

Parameters:

  • x (various shapes) [real/complex, inout/in,required] – in: x_in; out: mixed x

  • fx (various shapes) [real/complex, in] – x_out-x_in

  • alpha [real, in] – mixing parameter

interface  sf_optimize/adaptive_mix(x, fx, alpha, iter, z, fz)

Generic interface for the adaptive linear mixing, real (d_adaptive_mix) or complex (c_adaptive_mix): \(x \leftarrow x + \beta F_x\), with an element-wise mixing parameter \(\beta_j\), which starts from alpha, is increased by alpha at every iteration in which the residual \(F_{x,j}\) keeps its sign, up to 1, and is reset to alpha if the sign changes (oscillations). The previous residual is stored between calls: it is reset when iter is 1.

Parameters:

  • x (•) [real, inout/in,required] – in: x_in; out: mixed x

  • fx (size(x)) [real, in] – residual, x_out-x_in

  • alpha [real, in] – mixing parameter, initial and minimum speed

  • iter [integer, in] – iteration number, 1 resets the memory

  • z (•) [complex, inout/in,required] – z_in

  • fz (size(z)) [complex, in] – z_out-z_in

interface  sf_optimize/broyden_mix(x, fx, alpha, m, iter[, w0])

Generic interface for the modified Broyden mixing of Johnson (Phys. Rev. B 38, 12807, 1988), real (d_broyden_mix) or complex (c_broyden_mix): the mixed vector is built from the residual \(F_x = x_{out}-x_{in}\) and from the history of the last M iterations, stored between calls and reset when iter is 1. The first iteration is a linear mixing with alpha, as are all the iterations if M=0. The weight of the initial iteration is w0.

Parameters:

  • x (•) [real/complex, inout/in,required] – in: x_in; out: mixed x

  • fx (size(x)) [real/complex, in] – residual, x_out-x_in

  • alpha [real, in] – mixing parameter

  • m [integer, in] – number of previous iterations kept in memory, 0 is linear mixing

  • iter [integer, in] – iteration number, 1 resets the memory

Options:

w0 [real] – optional: weight of the initial iteration, default 0.01

interface  sf_optimize/fsolve(func, x, dfunc[, tol, info, icheck, maxfev])

Generic interface to solve a system of \(n\) nonlinear equations in \(n\) unknowns, \(f(x)=0\), with the hybrid Powell method of MINPACK. The Jacobian is computed numerically (hybrd, fsolve_hybrd_func and fsolve_hybrd_sub) or supplied by the user (hybrj, fsolve_hybrj_func and fsolve_hybrj_sub). The system is a function func(x) returning a vector with the size of x, or a subroutine func(x,f); the Jacobian dfunc(x) is a matrix of shape \([n,n]\). If icheck is true, the default, the program stops, writing the reason to a file, when the status returned by MINPACK is not a success.

Parameters:

  • func [external]

  • x (•) [real] – in: starting guess; out: solution, size n

  • dfunc [external]

Options:

  • tol [real] – optional: tolerance, default 1d-15

  • info [integer] – optional out: status of the MINPACK routine

  • icheck [logical] – optional: .true. (default) stops if the status is not a success

  • maxfev [integer] – optional: maximum number of function evaluations, default 200

subroutine  sf_optimize/brent(func, xmin[, brack, tol, niter])

This subroutine finds the minimum of a function of one variable, func, with the Brent method (golden section search and parabolic interpolation). The abscissa of the minimum is returned in xmin. The minimum must be bracketed: brack contains either a bracketing triplet (a,x,b), with a<x<b and func(x) lower than both func(a) and func(b), or two starting points, from which a bracketing triplet is searched downhill with bracket. If brack is absent the points 0 and 1 are used as starting points. The program stops if brack has a single element.

Parameters:

  • func [external] – function whose minimum is searched, real(8) function func(x) with real(8) x

  • xmin [real, inout] – abscissa of the minimum

Options:

  • brack (•) [real] – optional: 2 starting points or 3 bracketing points, default [0,1]

  • tol [real] – optional: fractional precision of the minimum, default 1d-9

  • niter [integer] – optional: maximum number of iterations, default 200

subroutine  sf_optimize/bracket(ax, bx, cx, fa, fb, fc, func)

Given a function func and two distinct initial points ax and bx, this subroutine searches in the downhill direction, defined by the function evaluated at the initial points, and returns three points ax, bx, cx which bracket a minimum of the function, i.e. bx is between ax and cx and func(bx) is lower than both func(ax) and func(cx). The values of the function at the three points are returned too.

Parameters:

  • ax [real, inout] – in: first initial point; out: first point of the bracketing triplet

  • bx [real, inout] – in: second initial point; out: middle point of the bracketing triplet

  • cx [real, out] – out: last point of the bracketing triplet

  • fa [real, out] – out: func(ax)

  • fb [real, out] – out: func(bx)

  • fc [real, out] – out: func(cx)

  • func [external] – function whose minimum is bracketed, real(8) func(x) with real(8) x

subroutine  sf_optimize/fmin(fn, start[, lambda, tol, conv_check, max_fun_calls, fun_calls, num_restart, ierr])

This subroutine minimizes a function of several variables with the Nelder-Mead downhill simplex algorithm (NELMIN), a simplex function minimisation procedure due to Nelder and Mead (1965), as implemented by O’Neill (Applied Statistics 20, 338, 1971). The convergence check, carried out every conv_check iterations, is on the variance of the function values at the vertices of the simplex; if the simplex is found not to be at a minimum it is restarted. The minimum value of the function is not returned.

Parameters:

  • fn [external] – function to be minimized, real(8) function fn(x) with real(8) x(:)

  • start (•) [real, in,required] – in: starting point; out: coordinates of the minimum

Options:

  • lambda (size(start)) [real] – optional: size and shape of the initial simplex, one step per variable, default 1

  • tol [real] – optional: terminating limit for the variance of the function values, > 0, default 1d-8

  • conv_check [integer] – optional: the convergence is checked every conv_check iterations, >= 1, default 10

  • max_fun_calls [integer] – optional: maximum number of function evaluations, default 500

  • fun_calls [integer] – optional out: number of function evaluations used

  • num_restart [integer] – optional out: number of restarts

  • ierr [integer] – optional out: 0 ok, 1 illegal tol or conv_check, 2 max_fun_calls exceeded

subroutine  sf_optimize/broyden1(ff, x[, icheck, maxits, tol, tol1, tolmin, stpmx, noexit])

This subroutine finds a zero of a system of \(n\) nonlinear equations in \(n\) unknowns, \(f(x)=0\), with the Broyden method and a backtracking line search (Numerical Recipes broydn): the Jacobian is computed numerically at the start, and then updated with the Broyden formula, restarting from a new numerical Jacobian when the update fails. On output x contains the solution. If the Jacobian turns singular the program stops, unless noexit is present and true.

Parameters:

  • ff [real] – system of equations, vector function of x

  • x (•) [real, inout] – in: starting guess; out: solution, size n

Options:

  • icheck [logical] – optional out: .true. if converged to a spurious minimum of \(|f|^2\)

  • maxits [integer] – optional: maximum number of iterations, default 200

  • tol [real] – optional: tolerance on the function values, default 1d-8

  • tol1 [real] – optional: tolerance at the starting point, default 0.01*tol

  • tolmin [real] – optional: spurious convergence tolerance, default 1d-7

  • stpmx [real] – optional: scaled maximum step of the line search, default 100

  • noexit [logical] – optional: .true. returns on singular Jacobian, set .false. if so

function  sf_optimize/brentq(func, a, b[, tol])

This function returns a zero of the function func in the interval [a,b], found with the Brent method. The function must change sign in the interval, i.e. the root must be bracketed: otherwise the program stops.

Parameters:

  • func [external] – function whose zero is searched, real(8) function func(x) with real(8) x

  • a [real, in] – lower end of the interval

  • b [real, in] – upper end of the interval

Options:

tol [real] – optional: tolerance, default epsilon(a)

Result:

fzero [real] – root of func in [a,b]

subroutine  sf_optimize/bisect(f, x1, x2[, eps, niter, flag])

This subroutine finds a zero of the function f in the interval [x1,x2] with the bisection method. On output x1 and x2 both contain the root. The flag flag is 0 if the function has the same sign at the two ends of the interval (no solution), the number of iterations if a root is found, and minus the number of iterations if the point found is a singularity rather than a root, i.e. \(|f|\ge 1\) there.

Parameters:

  • f [external] – function whose zero is searched, real(8) function f(x) with real(8) x

  • x1 [real] – in: lower end of the interval; out: root

  • x2 [real] – in: upper end of the interval; out: root

Options:

  • eps [real] – optional: width of the final interval, default 1d-9

  • niter [integer] – optional: maximum number of iterations, default 200

  • flag [integer] – optional out: >0 root found in flag iterations, 0 no solution, <0 singularity

subroutine  sf_optimize/fzero(f, b, c, iflag[, rguess, tol_rel, tol_abs])

This subroutine searches for a zero of the function f between the given values b and c, until the width of the interval collapses within the tolerance abs(b-c) <= 2*(rw*abs(b)+ae), with rw the relative and ae the absolute tolerance. The method is an efficient combination of bisection and the secant rule.

Control returns to the user in all cases, therefore iflag must be checked after each call. Its possible values are:

  • 1: b is within the requested tolerance of a zero. The interval (b,c) collapsed to the requested tolerance, the function changes sign in (b,c), and \(|f(x)|\) decreased as (b,c) collapsed.

  • 2: \(f(b)=0\). However, the interval (b,c) may not have collapsed to the requested tolerance.

  • 3: b may be near a singular point of \(f(x)\). The interval (b,c) collapsed to the requested tolerance and the function changes sign in (b,c), but \(|f(x)|\) increased as (b,c) collapsed, i.e. \(\max(|f(b_{in})|,|f(c_{in})|) < |f(b_{out})|\).

  • 4: no change in sign of \(f(x)\) was found although the interval (b,c) collapsed to the requested tolerance. The user must examine this case and decide whether b is near a local minimum of \(f(x)\), or near a zero of even multiplicity, or neither of these.

  • 5: too many (more than 500) function evaluations used.

Parameters:

  • f [external] – function whose zero is searched, real(8) function f(x) with real(8) x

  • b [real] – in: one end of the interval; out: best approximation of the zero

  • c [real] – in: other end of the interval; out: other end of the final interval

  • iflag [integer] – out: status code, see the list above

Options:

  • rguess [real] – optional: a better guess of the zero, default b; it is ignored if not interior to (b,c)

  • tol_rel [real] – optional: relative error rw of the stopping criterion; if smaller than machine precision it is set to it

  • tol_abs [real] – optional: absolute error ae of the stopping criterion; if (b,c) contains the origin it should be nonzero

subroutine  sf_optimize/newton(f, xinit[, eps, niter])

This subroutine finds a zero of the function f with the Newton method, starting from xinit, the derivative being evaluated numerically with finite differences. On output xinit contains the root.

Parameters:

  • f [external] – function whose zero is searched, real(8) function f(x) with real(8) x

  • xinit [real, inout] – in: starting point; out: root

Options:

  • eps [real] – optional: tolerance, default 1d-9

  • niter [integer] – optional: maximum number of iterations, default 200