.p>.CMCH8.DOC!BOUNDED_LEAST_SQUARES;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.
|
ibtype |
Action |
|
0 |
User will supply all the bounds. |
|
1 |
All variables are nonnegative |
|
2 |
All variables are nonpositive. |
|
3 |
User supplies only the bounds on 1st variable, all other variables will have the same bounds |
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 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
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 i-th component of the scaled
gradient at x is
calculated as
where g = ∇ F(x), s = xscale, and
Default: grad_tol = 
in double where ɛ is the machine precision
IMSL_STEP_TOL, float step_tol
(Input)
Scaled step tolerance. The i-th component of the scaled step
between two points x, and y, is computed as
where s = xscale.
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.
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 the Introduction, Passing Data to
User-Supplied Functions at the beginning of 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 the Introduction, Passing Data to
User-Supplied Functions at the beginning of 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 : Rn ® Rm, 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 + m I)-1 JTF
where m 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)|| £ e, li < xi < ui
g (xi) < 0, xi = ui
g (xi) >0, xi = li
where e 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).
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.
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 <math.h>
#define M 2
#define N 2
#define LDFJAC 2
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 <math.h>
#define M 2
#define N 2
#define LDFJAC 2
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
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |