# NAG FL Interfacee02gbf (glinc_​l1sol)

## ▸▿ Contents

Settings help

FL Name Style:

FL Specification Language:

## 1Purpose

e02gbf calculates an ${l}_{1}$ solution to an overdetermined system of linear equations, possibly subject to linear inequality constraints.

## 2Specification

Fortran Interface
 Subroutine e02gbf ( m, n, mpl, e, lde, f, x, mxs, k, el1n, indx, w, iw,
 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
#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)
The routine may be called by the names e02gbf or nagf_fit_glinc_l1sol.

## 3Description

Given a matrix $A$ with $m$ rows and $n$ columns $\left(m\ge n\right)$ and a vector $b$ with $m$ elements, the routine calculates an ${l}_{1}$ solution to the overdetermined system of equations
 $Ax=b.$
That is to say, it calculates a vector $x$, with $n$ elements, which minimizes the ${l}_{1}$-norm (the sum of the absolute values) of the residuals
 $r(x)=∑i=1m|ri|,$
where the residuals ${r}_{i}$ are given by
 $ri=bi-∑j=1naijxj, i=1,2,…,m.$
Here ${a}_{ij}$ is the element in row $i$ and column $j$ of $A$, ${b}_{i}$ is the $i$th element of $b$ and ${x}_{j}$ the $j$th element of $x$.
If, in addition, a matrix $C$ with $l$ rows and $n$ columns and a vector $d$ with $l$ elements, are given, the vector $x$ computed by the routine is such as to minimize the ${l}_{1}$-norm $r\left(x\right)$ subject to the set of inequality constraints $Cx\ge d$.
The matrices $A$ and $C$ need not be of full rank.
Typically in applications to data fitting, data consisting of $m$ points with coordinates $\left({t}_{i},{y}_{i}\right)$ is to be approximated by a linear combination of known functions ${\varphi }_{i}\left(t\right)$,
 $α1ϕ1(t)+α2ϕ2(t)+⋯+αnϕn(t),$
in the ${l}_{1}$-norm, possibly subject to linear inequality constraints on the coefficients ${\alpha }_{j}$ of the form $C\alpha \ge d$ where $\alpha$ is the vector of the ${\alpha }_{j}$ and $C$ and $d$ are as in the previous paragraph. This is equivalent to finding an ${l}_{1}$ solution to the overdetermined system of equations
 $∑j=1nϕj(ti)αj=yi, i=1,2,…,m,$
subject to $C\alpha \ge d$.
Thus if, for each value of $i$ and $j$, the element ${a}_{ij}$ of the matrix $A$ above is set equal to the value of ${\varphi }_{j}\left({t}_{i}\right)$ and ${b}_{i}$ is equal to ${y}_{i}$ and $C$ and $d$ are also supplied to the routine, the solution vector $x$ will contain the required values of the ${\alpha }_{j}$. Note that the independent variable $t$ above can, instead, be a vector of several independent variables (this includes the case where each of ${\varphi }_{i}$ 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
 $g(x) = γ r(x) - ∑ i=1 l min(0, ciT x-di ) ,$
where $\gamma$ is a penalty parameter, ${c}_{i}^{\mathrm{T}}$ is the $i$th row of the matrix $C$, and ${d}_{i}$ is the $i$th element of the vector $d$. 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.

## 4References

Bartels R H, Conn A R and Charalambous C (1976) Minimisation techniques for piecewise Differentiable functions – the ${l}_{\infty }$ 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 ${l}_{1}$ 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 ${l}_{1}$ 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

## 5Arguments

1: $\mathbf{m}$Integer Input
On entry: the number of equations in the overdetermined system, $m$ (i.e., the number of rows of the matrix $A$).
Constraint: ${\mathbf{m}}\ge {\mathbf{n}}$.
2: $\mathbf{n}$Integer Input
On entry: the number of unknowns, $n$ (the number of columns of the matrix $A$).
Constraint: ${\mathbf{n}}\ge 2$.
3: $\mathbf{mpl}$Integer Input
On entry: $m+l$, where $l$ is the number of constraints (which may be zero).
Constraint: ${\mathbf{mpl}}\ge {\mathbf{m}}$.
4: $\mathbf{e}\left({\mathbf{lde}},{\mathbf{mpl}}\right)$Real (Kind=nag_wp) array Input/Output
On entry: the equation and constraint matrices stored in the following manner.
The first $m$ columns contain the $m$ rows of the matrix $A$; element ${\mathbf{e}}\left(\mathit{i},\mathit{j}\right)$ specifying the element ${a}_{\mathit{j}\mathit{i}}$ in the $\mathit{j}$th row and $\mathit{i}$th column of $A$ (the coefficient of the $\mathit{i}$th unknown in the $\mathit{j}$th equation), for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,m$. The next $l$ columns contain the $l$ rows of the constraint matrix $C$; element ${\mathbf{e}}\left(\mathit{i},\mathit{j}+m\right)$ containing the element ${c}_{\mathit{j}\mathit{i}}$ in the $\mathit{j}$th row and $\mathit{i}$th column of $C$ (the coefficient of the $\mathit{i}$th unknown in the $\mathit{j}$th constraint), for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,l$.
On exit: unchanged, except possibly to the extent of a small multiple of the machine precision. (See Section 9.)
5: $\mathbf{lde}$Integer Input
On entry: the first dimension of the array e as declared in the (sub)program from which e02gbf is called.
Constraint: ${\mathbf{lde}}\ge {\mathbf{n}}$.
6: $\mathbf{f}\left({\mathbf{mpl}}\right)$Real (Kind=nag_wp) array Input
On entry: ${\mathbf{f}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,m$, must contain ${b}_{\mathit{i}}$ (the $\mathit{i}$th element of the right-hand side vector of the overdetermined system of equations) and ${\mathbf{f}}\left(m+\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,l$, must contain ${d}_{i}$ (the $i$th element of the right-hand side vector of the constraints), where $l$ is the number of constraints.
7: $\mathbf{x}\left({\mathbf{n}}\right)$Real (Kind=nag_wp) array Input/Output
On entry: ${\mathbf{x}}\left(\mathit{i}\right)$ must contain an estimate of the $\mathit{i}$th unknown, for $\mathit{i}=1,2,\dots ,n$. If no better initial estimate for ${\mathbf{x}}\left(i\right)$ is available, set ${\mathbf{x}}\left(i\right)=0.0$.
On exit: the latest estimate of the $\mathit{i}$th unknown, for $\mathit{i}=1,2,\dots ,n$. If ${\mathbf{ifail}}={\mathbf{0}}$ on exit, these are the solution values.
8: $\mathbf{mxs}$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 $n$. If, on entry, mxs is zero or negative, the value returned by x02bbf is used.
9: $\mathbf{monit}$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
 Subroutine monit ( n, x, k, el1n)
 Integer, Intent (In) :: n, niter, k Real (Kind=nag_wp), Intent (In) :: x(n), el1n
 void monit (const Integer *n, const double x[], const Integer *niter, const Integer *k, const double *el1n)
1: $\mathbf{n}$Integer Input
On entry: the number $n$ of unknowns (the number of columns of the matrix $A$).
2: $\mathbf{x}\left({\mathbf{n}}\right)$Real (Kind=nag_wp) array Input
On entry: the latest estimate of the unknowns.
3: $\mathbf{niter}$Integer Input
On entry: the number of iterations so far carried out.
4: $\mathbf{k}$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: $\mathbf{el1n}$Real (Kind=nag_wp) Input
On entry: the ${l}_{1}$-norm of the current residuals of the overdetermined 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: $\mathbf{iprint}$Integer Input
On entry: the frequency of iteration print out.
${\mathbf{iprint}}>0$
monit is called every iprint iterations and at the solution.
${\mathbf{iprint}}=0$
Information is printed out at the solution only. Otherwise monit is not called (but a dummy routine must still be provided).
11: $\mathbf{k}$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: $\mathbf{el1n}$Real (Kind=nag_wp) Output
On exit: the ${l}_{1}$-norm (sum of absolute values) of the equation residuals.
13: $\mathbf{indx}\left({\mathbf{mpl}}\right)$Integer array Output
On exit: specifies which columns of e relate to the inactive equations and constraints. ${\mathbf{indx}}\left(1\right)$ up to ${\mathbf{indx}}\left({\mathbf{k}}\right)$ number the active columns and ${\mathbf{indx}}\left({\mathbf{k}}+1\right)$ up to ${\mathbf{indx}}\left({\mathbf{mpl}}\right)$ number the inactive columns.
14: $\mathbf{w}\left({\mathbf{iw}}\right)$Real (Kind=nag_wp) array Workspace
15: $\mathbf{iw}$Integer Input
On entry: the dimension of the array w as declared in the (sub)program from which e02gbf is called.
Constraint: ${\mathbf{iw}}\ge 3×{\mathbf{mpl}}+5×{\mathbf{n}}+{{\mathbf{n}}}^{2}+\left({\mathbf{n}}+1\right)×\left({\mathbf{n}}+2\right)/2$.
16: $\mathbf{ifail}$Integer Input/Output
On entry: ifail must be set to $0$, $-1$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $-1$ means that an error message is printed while a value of $1$ means that it is not.
If halting is not appropriate, the value $-1$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $0$ is recommended. When the value $-\mathbf{1}$ or $\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit: ${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

## 6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-1$, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
${\mathbf{ifail}}=1$
The constraints cannot all be satisfied simultaneously: they are not compatible with one another. Hence no solution is possible.
${\mathbf{ifail}}=2$
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.
${\mathbf{ifail}}=3$
The routine has failed because of numerical difficulties; the problem is too ill-conditioned. Consider rescaling the unknowns.
${\mathbf{ifail}}=4$
Elements $1$ 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 $A$ or $C$.
On entry, iw is too small. ${\mathbf{iw}}=⟨\mathit{\text{value}}⟩$. Minimum possible dimension: $⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{lde}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{lde}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{m}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{m}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{mpl}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{m}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{mpl}}\ge {\mathbf{m}}$.
On entry, ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{n}}\ge 2$.
${\mathbf{ifail}}=-99$
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{ifail}}=-399$
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.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

## 7Accuracy

The method is stable.

## 8Parallelism 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 $m$ and $n$ 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 $n$ and the total time taken is approximately proportional to $m{n}^{2}$.
Linear dependencies among the rows or columns of $A$ and $C$ are not necessarily a problem to the algorithm. Solutions can be obtained from rank-deficient $A$ and $C$. 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.

## 10Example

Suppose we wish to approximate in $\left[0,1\right]$ a set of data by a curve of the form
 $y=ax3+bx2+cx+d$
which has non-negative slope at the data points. Given points $\left({t}_{i},{y}_{i}\right)$ we may form the equations
 $yi=ati3+bti2+cti+d$
for $\mathit{i}=1,2,\dots ,6$, for the $6$ data points. The requirement of a non-negative slope at the data points demands
 $3ati2+2bti+c≥0$
for each ${t}_{i}$ 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.1Program Text

Program Text (e02gbfe.f90)

### 10.2Program Data

Program Data (e02gbfe.d)

### 10.3Program Results

Program Results (e02gbfe.r)