NAG FL Interface
e02gbf (glinc_l1sol)
1
Purpose
e02gbf calculates an solution to an over-determined system of linear equations, possibly subject to linear inequality constraints.
2
Specification
Fortran Interface
Subroutine e02gbf ( |
m, n, mpl, e, lde, f, x, mxs, monit, iprint, k, el1n, indx, w, iw, ifail) |
Integer, Intent (In) |
:: |
m, n, mpl, lde, mxs, iprint, iw |
Integer, Intent (Inout) |
:: |
ifail |
Integer, Intent (Out) |
:: |
k, indx(mpl) |
Real (Kind=nag_wp), Intent (In) |
:: |
f(mpl) |
Real (Kind=nag_wp), Intent (Inout) |
:: |
e(lde,mpl), x(n) |
Real (Kind=nag_wp), Intent (Out) |
:: |
el1n, w(iw) |
External |
:: |
monit |
|
C Header Interface
#include <nag.h>
void |
e02gbf_ (const Integer *m, const Integer *n, const Integer *mpl, double e[], const Integer *lde, const double f[], double x[], const Integer *mxs, void (NAG_CALL *monit)(const Integer *n, const double x[], const Integer *niter, const Integer *k, const double *el1n), const Integer *iprint, Integer *k, double *el1n, Integer indx[], double w[], const Integer *iw, Integer *ifail) |
|
C++ Header Interface
#include <nag.h> extern "C" {
void |
e02gbf_ (const Integer &m, const Integer &n, const Integer &mpl, double e[], const Integer &lde, const double f[], double x[], const Integer &mxs, void (NAG_CALL *monit)(const Integer &n, const double x[], const Integer &niter, const Integer &k, const double &el1n), const Integer &iprint, Integer &k, double &el1n, Integer indx[], double w[], const Integer &iw, Integer &ifail) |
}
|
The routine may be called by the names e02gbf or nagf_fit_glinc_l1sol.
3
Description
Given a matrix
with
rows and
columns
and a vector
with
elements, the routine calculates an
solution to the over-determined system of equations
That is to say, it calculates a vector
, with
elements, which minimizes the
-norm (the sum of the absolute values) of the residuals
where the residuals
are given by
Here
is the element in row
and column
of
,
is the
th element of
and
the
th element of
.
If, in addition, a matrix with rows and columns and a vector with elements, are given, the vector computed by the routine is such as to minimize the -norm subject to the set of inequality constraints .
The matrices and need not be of full rank.
Typically in applications to data fitting, data consisting of
points with coordinates
is to be approximated by a linear combination of known functions
,
in the
-norm, possibly subject to linear inequality constraints on the coefficients
of the form
where
is the vector of the
and
and
are as in the previous paragraph. This is equivalent to finding an
solution to the over-determined system of equations
subject to
.
Thus if, for each value of and , the element of the matrix above is set equal to the value of and is equal to and and are also supplied to the routine, the solution vector will contain the required values of the . Note that the independent variable above can, instead, be a vector of several independent variables (this includes the case where each of is a function of a different variable, or set of variables).
The algorithm follows the Conn–Pietrzykowski approach (see
Bartels et al. (1978) and
Conn and Pietrzykowski (1977)), which is via an exact penalty function
where
is a penalty parameter,
is the
th row of the matrix
, and
is the
th element of the vector
. It proceeds in a step-by-step manner much like the simplex method for linear programming but does not move from vertex to vertex and does not require the problem to be cast in a form containing only non-negative unknowns. It uses stable procedures to update an orthogonal factorization of the current set of active equations and constraints.
4
References
Bartels R H, Conn A R and Charalambous C (1976) Minimisation techniques for piecewise Differentiable functions – the solution to an overdetermined linear system Technical Report No. 247, CORR 76/30 Mathematical Sciences Department, The John Hopkins University
Bartels R H, Conn A R and Sinclair J W (1976) A Fortran program for solving overdetermined systems of linear equations in the Sense Technical Report No. 236, CORR 76/7 Mathematical Sciences Department, The John Hopkins University
Bartels R H, Conn A R and Sinclair J W (1978) Minimisation techniques for piecewise differentiable functions – the solution to an overdetermined linear system SIAM J. Numer. Anal. 15 224–241
Conn A R and Pietrzykowski T (1977) A penalty-function method converging directly to a constrained optimum SIAM J. Numer. Anal. 14 348–375
5
Arguments
-
1:
– Integer
Input
-
On entry: the number of equations in the over-determined system, (i.e., the number of rows of the matrix ).
Constraint:
.
-
2:
– Integer
Input
-
On entry: the number of unknowns, (the number of columns of the matrix ).
Constraint:
.
-
3:
– Integer
Input
-
On entry: , where is the number of constraints (which may be zero).
Constraint:
.
-
4:
– Real (Kind=nag_wp) array
Input/Output
-
On entry: the equation and constraint matrices stored in the following manner.
The first columns contain the rows of the matrix ; element
specifying the element in the th row and th column of (the coefficient of the th unknown in the th equation), for and . The next columns contain the rows of the constraint matrix ; element
containing the element in the th row and th column of (the coefficient of the th unknown in the th constraint), for and .
On exit: unchanged, except possibly to the extent of a small multiple of the
machine precision. (See
Section 9.)
-
5:
– Integer
Input
-
On entry: the first dimension of the array
e as declared in the (sub)program from which
e02gbf is called.
Constraint:
.
-
6:
– Real (Kind=nag_wp) array
Input
-
On entry: , for , must contain
(the th element of the right-hand side vector of the over-determined system of equations) and , for , must contain (the th element of the right-hand side vector of the constraints), where is the number of constraints.
-
7:
– Real (Kind=nag_wp) array
Input/Output
-
On entry: must contain an estimate of the th unknown, for . If no better initial estimate for is available, set .
On exit: the latest estimate of the
th unknown, for . If on exit, these are the solution values.
-
8:
– Integer
Input
-
On entry: the maximum number of steps to be allowed for the solution of the unconstrained problem. Typically this may be a modest multiple of
. If, on entry,
mxs is zero or negative, the value returned by
x02bbf is used.
-
9:
– Subroutine, supplied by the user.
External Procedure
-
monit can be used to print out the current values of any selection of its arguments. The frequency with which
monit is called in
e02gbf is controlled by
iprint.
The specification of
monit is:
Fortran Interface
Integer, Intent (In) |
:: |
n, niter, k |
Real (Kind=nag_wp), Intent (In) |
:: |
x(n), el1n |
|
C Header Interface
void |
monit_ (const Integer *n, const double x[], const Integer *niter, const Integer *k, const double *el1n) |
|
C++ Header Interface
#include <nag.h> extern "C" {
void |
monit_ (const Integer &n, const double x[], const Integer &niter, const Integer &k, const double &el1n) |
}
|
-
1:
– Integer
Input
-
On entry: the number of unknowns (the number of columns of the matrix ).
-
2:
– Real (Kind=nag_wp) array
Input
-
On entry: the latest estimate of the unknowns.
-
3:
– Integer
Input
-
On entry: the number of iterations so far carried out.
-
4:
– Integer
Input
-
On entry: the total number of equations and constraints which are currently active (i.e., the number of equations with zero residuals plus the number of constraints which are satisfied as equations).
-
5:
– Real (Kind=nag_wp)
Input
-
On entry: the -norm of the current residuals of the over-determined system of equations.
monit must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
e02gbf is called. Arguments denoted as
Input must
not be changed by this procedure.
-
10:
– Integer
Input
-
On entry: the frequency of iteration print out.
- monit is called every iprint iterations and at the solution.
- Information is printed out at the solution only. Otherwise monit is not called (but a dummy routine must still be provided).
-
11:
– Integer
Output
-
On exit: the total number of equations and constraints which are then active (i.e., the number of equations with zero residuals plus the number of constraints which are satisfied as equalities).
-
12:
– Real (Kind=nag_wp)
Output
-
On exit: the -norm (sum of absolute values) of the equation residuals.
-
13:
– Integer array
Output
-
On exit: specifies which columns of
e relate to the inactive equations and constraints.
up to
number the active columns and
up to
number the inactive columns.
-
14:
– Real (Kind=nag_wp) array
Workspace
-
15:
– Integer
Input
-
On entry: the dimension of the array
w as declared in the (sub)program from which
e02gbf is called.
Constraint:
.
-
16:
– Integer
Input/Output
-
On entry:
ifail must be set to
,
. If you are unfamiliar with this argument you should refer to
Section 4 in the Introduction to the NAG Library FL Interface for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
.
When the value is used it is essential to test the value of ifail on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
-
The constraints cannot all be satisfied simultaneously: they are not compatible with one another. Hence no solution is possible.
-
The limit imposed by
mxs has been reached without finding a solution. Consider restarting from the current point by simply calling
e02gbf again without changing the arguments.
-
The routine has failed because of numerical difficulties; the problem is too ill-conditioned. Consider rescaling the unknowns.
-
Elements
to
m of one of the first
mpl columns of the array
e are all zero – this corresponds to a zero row in either of the matrices
or
.
On entry,
iw is too small.
. Minimum possible dimension:
.
On entry, and .
Constraint: .
On entry, and .
Constraint: .
On entry, and .
Constraint: .
On entry, .
Constraint: .
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 7 in the Introduction to the NAG Library FL Interface for further information.
Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library FL Interface for further information.
Dynamic memory allocation failed.
See
Section 9 in the Introduction to the NAG Library FL Interface for further information.
7
Accuracy
The method is stable.
8
Parallelism and Performance
e02gbf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
e02gbf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementation-specific information.
The effect of and on the time and on the number of iterations varies from problem to problem, but typically the number of iterations is a small multiple of and the total time taken is approximately proportional to .
Linear dependencies among the rows or columns of
and
are not necessarily a problem to the algorithm. Solutions can be obtained from rank-deficient
and
. However, the algorithm requires that at every step the currently active columns of
e form a linearly independent set. If this is not the case at any step, small, random perturbations of the order of rounding error are added to the appropriate columns of
e. Normally this perturbation process will not affect the solution significantly. It does mean, however, that results may not be exactly reproducible.
10
Example
Suppose we wish to approximate in
a set of data by a curve of the form
which has non-negative slope at the data points. Given points
we may form the equations
for
, for the
data points. The requirement of a non-negative slope at the data points demands
for each
and these form the constraints.
(Note that, for fitting with polynomials, it would usually be advisable to work with the polynomial expressed in Chebyshev series form (see the
E02 Chapter Introduction). The power series form is used here for simplicity of exposition.)
10.1
Program Text
10.2
Program Data
10.3
Program Results