bounded_least_squares

Solves a nonlinear least-squares problem subject to bounds on the variables using a modified Levenberg-Marquardt algorithm.

Synopsis

#include <imsl.h>

float *imsl_f_bounded_least_squares (void fcn(), int m, int n, int ibtype, float xlb[], float xub[], ..., 0)

The type double function is imsl_d_bounded_least_squares.

Required Arguments

void fcn (int m, int n, float x[], float f[]) (Input/Output)
User-supplied function to evaluate the function that defines the least-squares problem where x is a vector of length n at which point the function is evaluated, and f is a vector of length m containing the function values at point x.

int m (Input)
Number of functions.

int n (Input)
Number of variables where n ≤ m.

int ibtype (Input)
Scalar indicating the types of bounds on the variables.

float xlb[] (Input, Output, or Input/Output)
Array with n components containing the lower bounds on the variables. (Input, if ibtype = 0; output, if ibtype = 1 or 2; Input/Output, if ibtype = 3)

If there is no lower bound on a variable, then the corresponding xlb value should be set to −106.

float xub[] (Input, Output, or Input/Output)
Array with n components containing the upper bounds on the variables. (Input, if ibtype = 0; output, if ibtype 1 or 2; Input/Output, if ibtype = 3)

If there is no upper bound on a variable, then the corresponding xub value should be set to 106.

Return Value

A pointer to the solution x of the nonlinear least-squares problem. To release this space, use imsl_free. If no solution can be computed, then NULL is returned.

Synopsis with Optional Arguments

#include <imsl.h>

float *imsl_f_bounded_least_squares (void fcn(), int m, int n, int ibtype, float xlb[], float xub[],

IMSL_XGUESS, float xguess[],

IMSL_JACOBIAN, void jacobian(),

IMSL_XSCALE, float xscale[],

IMSL_FSCALE, float fscale[],

IMSL_GRAD_TOL, float grad_tol,

IMSL_STEP_TOL, float step_tol,

IMSL_REL_FCN_TOL, float rfcn_tol,

IMSL_ABS_FCN_TOL, float afcn_tol,

IMSL_MAX_STEP, float max_step,

IMSL_INIT_TRUST_REGION, float trust_region,

IMSL_GOOD_DIGIT, int ndigit,

IMSL_MAX_ITN, int max_itn,

IMSL_MAX_FCN, int max_fcn,

IMSL_MAX_JACOBIAN, int max_jacobian,

IMSL_INTERN_SCALE,

IMSL_RETURN_USER, float x[],

IMSL_FVEC, float **fvec,

IMSL_FVEC_USER, float fvec[],

IMSL_FJAC, float **fjac,

IMSL_FJAC_USER, float fjac[],

IMSL_FJAC_COL_DIM, int fjac_col_dim,

IMSL_FCN_W_DATA, void fcn(), void *data,

IMSL_JACOBIAN_W_DATA, void jacobian(), void *data,

0)

Optional Arguments

IMSL_XGUESS, float xguess[] (Input)
Array with n components containing an initial guess.
Default: xguess = 0

IMSL_JACOBIAN, void jacobian (int m, int n, float x[], float fjac[], int fjac_col_dim) (Input)
User-supplied function to compute the Jacobian where x is a vector of length n at which point the Jacobian is evaluated, fjac is the computed m × n Jacobian at the point x, and fjac_col_dim is the column dimension of fjac. Note that each partial derivative ∂fi /∂xj should be returned in fjac[(i‑1)*fjac_col_dim+j-1].

IMSL_XSCALE, float xscale[] (Input)
Array with n components containing the scaling vector for the variables. Argument xscale is used mainly in scaling the gradient and the distance between two points. See keywords IMSL_GRAD_TOL and IMSL_STEP_TOL for more details.
Default: xscale[] = 1

IMSL_FSCALE, float fscale[] (Input)
Array with m components containing the diagonal scaling matrix for the functions. The i-th component of fscale is a positive scalar specifying the reciprocal magnitude of the i-th component function of the problem.
Default: fscale[] = 1

IMSL_GRAD_TOL, float grad_tol (Input)
Scaled gradient tolerance.

The second bounded_least_squares stopping criterion occurs when

≤ grad_tol

where gsi(x) is component i of the scaled gradient of F at x, defined as:

,

and where J is the Jacobian matrix for m-component function vector F(x) with n‑component argument x with Jji = ∇ifj(x), si = xscale[i-1], and fsi = fscale[i-1].

Default: grad_tol =  in double where ɛ is the machine precision.

IMSL_STEP_TOL, float step_tol (Input)
Scaled step tolerance.

The third bounded_least_squares stopping criterion occurs when

 ≤ step_tol,

where Δxi denotes the i-th component of the iteration step Δx,xi is the i-th component of the current iterate x, and si = xscale[i-1].

Default: step_tol = ɛ2/3, where ɛ is the machine precision

IMSL_REL_FCN_TOL, float rfcn_tol (Input)
Relative function tolerance.

Default: rfcn_tol = max(10-10, ɛ2/3), max(10-20, ɛ2/3) in double, where ɛ is the machine precision

IMSL_ABS_FCN_TOL, float afcn_tol (Input)
Absolute function tolerance.

The first bounded_least_squares stopping criterion occurs when objective function

≤ afcn_tol.

Default: afcn_tol = max(10-20, ɛ2), max(10-40, ɛ2) in double, where ɛ is the machine precision

IMSL_MAX_STEP, float max_step (Input)
Maximum allowable step size.
Default: max_step = 1000 max(ɛ1, ɛ2), where

for s = xscale and t = xguess.

IMSL_INIT_TRUST_REGION, float trust_region (Input)
Size of initial trust region radius. The default is based on the initial scaled Cauchy step.

IMSL_GOOD_DIGIT, int ndigit (Input)
Number of good digits in the function.
Default: machine dependent

IMSL_MAX_ITN, int max_itn (Input)
Maximum number of iterations.
Default: max_itn = 100

IMSL_MAX_FCN, int max_fcn (Input)
Maximum number of function evaluations.
Default: max_fcn = 400

IMSL_MAX_JACOBIAN, int max_jacobian (Input)
Maximum number of Jacobian evaluations.
Default: max_jacobian = 400

IMSL_INTERN_SCALE
Internal variable scaling option. With this option, the values for xscale are set internally.

IMSL_RETURN_USER, float x[] (Output)
Array with n components containing the computed solution.

IMSL_FVEC, float **fvec (Output)
The address of a pointer to a real array of length m containing the residuals at the approximate solution. On return, the necessary space is allocated by imsl_f_bounded_least_squares. Typically, float *fvec is declared, and &fvec is used as an argument.

IMSL_FVEC_USER, float fvec[] (Output)
A user-allocated array of size m containing the residuals at the approximate solution.

IMSL_FJAC, float **fjac (Output)
The address of a pointer to an array of size m × n containing the Jacobian at the approximate solution. On return, the necessary space is allocated by imsl_f_bounded_least_squares. Typically, float *fjac is declared, and &fjac is used as an argument.

IMSL_FJAC_USER, float fjac[] (Output)
A user-allocated array of size m × n containing the Jacobian at the approximate solution.

IMSL_FJAC_COL_DIM, int fjac_col_dim (Input)
The column dimension of fjac.
Default: fjac_col_dim = n

IMSL_FCN_W_DATA, void fcn (int m, int n, float x[], float f[], void *data), void *data, (Input)
User-supplied function to evaluate the function that defines the least-squares problem, which also accepts a pointer to data that is supplied by the user. data is a pointer to the data to be passed to the user-supplied function. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.

IMSL_JACOBIAN_W_DATA, void jacobian (int m, int n, float x[], float fjac[], int fjac_col_dim, void *data), void *data, (Input)
User-supplied function to compute the Jacobian, which also accepts a pointer to data that is supplied by the user. data is a pointer to the data to be passed to the user-supplied function. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.

Description

The function imsl_f_bounded_least_squares uses a modified Levenberg-Marquardt method and an active set strategy to solve nonlinear least-squares problems subject to simple bounds on the variables. The problem is stated as follows:

subject to l ≤x ≤u

where m ≥n, F : ℜn → ℜm, and fi(x) is the i-th component function of F(x). From a given starting point, an active set IA, which contains the indices of the variables at their bounds, is built. A variable is called a “free variable” if it is not in the active set. The routine then computes the search direction for the free variables according to the formula

d =-(JTJ +μ I)-1 JTF

where μ is the Levenberg-Marquardt parameter, F =F(x), and J is the Jacobian with respect to the free variables. The search direction for the variables in IA is set to zero. The trust region approach discussed by Dennis and Schnabel (1983) is used to find the new point. Finally, the optimality conditions are checked. The conditions are

∥g (xi)∥ ≤ɛ, li < xi < ui

g (xi) < 0, xi =ui

g (xi) >0, xi =li

where ɛ is a gradient tolerance. This process is repeated until the optimality criterion is achieved.

The active set is changed only when a free variable hits its bounds during an iteration or the optimality condition is met for the free variables but not for all variables in IA, the active set. In the latter case, a variable that violates the optimality condition will be dropped out of IA. For more detail on the Levenberg-Marquardt method, see Levenberg (1944) or Marquardt (1963). For more detail on the active set strategy, see Gill and Murray (1976).

The first stopping criterion for imsl_f_bounded_least_squares occurs when the norm of the function is less than the absolute function tolerance. The second stopping criterion occurs when the norm of the scaled gradient is less than the scaled gradient tolerance. The third stopping criterion occurs when the scaled distance between the last two steps is less than the step tolerance. See options IMSL_ABS_FCN_TOL, IMSL_GRAD_TOL, and IMSL_STEP_TOL for details.

Since a finite-difference method is used to estimate the Jacobian for some single-precision calculations, an inaccurate estimate of the Jacobian may cause the algorithm to terminate at a noncritical point. In such cases, high-precision arithmetic is recommended. Also, whenever the exact Jacobian can be easily provided, the option IMSL_JACOBIAN should be used.

On some platforms, imsl_f_bounded_least_squares can evaluate the user-supplied functions fcn and jacobian in parallel. This is done only if the function imsl_omp_options is called to flag user-defined functions as thread-safe. A function is thread-safe if there are no dependencies between calls. Such dependencies are usually the result of writing to global or static variables.

Examples

Example 1

In this example, the nonlinear least-squares problem

where

is solved with an initial guess (-1.2, 1.0).

#include <imsl.h>

#include <stdio.h>

#include <math.h>

#define M 2

#define N 2

#define LDFJAC 2

int main()

{

void rosbck(int, int, float *, float *);

int ibtype = 0;

float xlb[N] = {-2.0, -1.0};

float xub[N] = {0.5, 2.0};

float *x;

x = imsl_f_bounded_least_squares (rosbck, M, N, ibtype, xlb,

xub, 0);

printf("x[0] = %f\n", x[0]);

printf("x[1] = %f\n", x[1]);

}

void rosbck (int m, int n, float *x, float *f)

{

f[0] = 10.0*(x[1] - x[0]*x[0]);

f[1] = 1.0 - x[0];

}

Output

x[0] = 0.500000
x[1] = 0.250000

Example 2

This example solves the nonlinear least-squares problem

where

This time, an initial guess (-1.2, 1.0) is supplied, as well as the analytic Jacobian. The residual at the approximate solution is returned.

#include <imsl.h>

#include <stdio.h>

#include <math.h>

#define M 2

#define N 2

#define LDFJAC 2

int main()

{

void rosbck(int, int, float *, float *);

void jacobian(int, int, float *, float *, int);

int ibtype = 0;

float xlb[N] = {-2.0, -1.0};

float xub[N] = {0.5, 2.0};

float xguess[N] = {-1.2, 1.0};

float *fvec;

float *x;

x = imsl_f_bounded_least_squares (rosbck, M, N, ibtype, xlb, xub,

IMSL_JACOBIAN, jacobian,

IMSL_XGUESS, xguess,

IMSL_FVEC, &fvec,

0);

printf("x[0] = %f\n", x[0]);

printf("x[1] = %f\n\n", x[1]);

printf("fvec[0] = %f\n", fvec[0]);

printf("fvec[1] = %f\n\n", fvec[1]);

}

void rosbck (int m, int n, float *x, float *f)

{

f[0] = 10.0*(x[1] - x[0]*x[0]);

f[1] = 1.0 - x[0];

}

void jacobian (int m, int n, float *x, float *fjac, int fjac_col_dim)

{

fjac[0] = -20.0*x[0];

fjac[1] = 10.0;

fjac[2] = -1.0;

fjac[3] = 0.0;

}

Output

x[0] = 0.500000

x[1] = 0.250000

fvec[0] = 0.000000

fvec[1] = 0.500000

Fatal Errors