PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_fit_glinc_l1sol (e02gb)
Purpose
nag_fit_glinc_l1sol (e02gb) calculates an ${l}_{1}$ solution to an overdetermined system of linear equations, possibly subject to linear inequality constraints.
Syntax
[
e,
x,
k,
el1n,
indx,
ifail] = e02gb(
m,
e,
f,
x,
mxs,
monit,
iprint, 'n',
n, 'mpl',
mpl)
[
e,
x,
k,
el1n,
indx,
ifail] = nag_fit_glinc_l1sol(
m,
e,
f,
x,
mxs,
monit,
iprint, 'n',
n, 'mpl',
mpl)
Description
Given a matrix
$A$ with
$m$ rows and
$n$ columns
$\left(m\ge n\right)$ and a vector
$b$ with
$m$ elements, the function calculates an
${l}_{1}$ solution to the overdetermined system of equations
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
where the residuals
${r}_{i}$ are given by
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 function 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)$,
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
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 function, 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
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 stepbystep 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 nonnegative unknowns. It uses stable procedures to update an orthogonal factorization of the current set of active equations and constraints.
References
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 penaltyfunction method converging directly to a constrained optimum SIAM J. Numer. Anal. 14 348–375
Parameters
Compulsory Input Parameters
 1:
$\mathrm{m}$ – int64int32nag_int scalar

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:
$\mathrm{e}\left(\mathit{lde},{\mathbf{mpl}}\right)$ – double array

lde, the first dimension of the array, must satisfy the constraint
$\mathit{lde}\ge {\mathbf{n}}$.
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$.
 3:
$\mathrm{f}\left({\mathbf{mpl}}\right)$ – double array

${\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 righthand 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 righthand side vector of the constraints), where $l$ is the number of constraints.
 4:
$\mathrm{x}\left({\mathbf{n}}\right)$ – double array

${\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$.
 5:
$\mathrm{mxs}$ – int64int32nag_int scalar

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
nag_machine_integer_max (x02bb) is used.
 6:
$\mathrm{monit}$ – function handle or string containing name of mfile

monit can be used to print out the current values of any selection of its arguments. The frequency with which
monit is called in
nag_fit_glinc_l1sol (e02gb) is controlled by
iprint.
monit(n, x, niter, k, el1n)
Input Parameters
 1:
$\mathrm{n}$ – int64int32nag_int scalar

The number $n$ of unknowns (the number of columns of the matrix $A$).
 2:
$\mathrm{x}\left({\mathbf{n}}\right)$ – double array

The latest estimate of the unknowns.
 3:
$\mathrm{niter}$ – int64int32nag_int scalar

The number of iterations so far carried out.
 4:
$\mathrm{k}$ – int64int32nag_int scalar

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:
$\mathrm{el1n}$ – double scalar

The ${l}_{1}$norm of the current residuals of the overdetermined system of equations.
 7:
$\mathrm{iprint}$ – int64int32nag_int scalar

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 function must still be provided).
Optional Input Parameters
 1:
$\mathrm{n}$ – int64int32nag_int scalar

Default:
the dimension of the array
x and the first dimension of the array
e. (An error is raised if these dimensions are not equal.)
The number of unknowns, $n$ (the number of columns of the matrix $A$).
Constraint:
${\mathbf{n}}\ge 2$.
 2:
$\mathrm{mpl}$ – int64int32nag_int scalar

Default:
the dimension of the array
f and the second dimension of the array
e. (An error is raised if these dimensions are not equal.)
$m+l$, where $l$ is the number of constraints (which may be zero).
Constraint:
${\mathbf{mpl}}\ge {\mathbf{m}}$.
Output Parameters
 1:
$\mathrm{e}\left(\mathit{lde},{\mathbf{mpl}}\right)$ – double array

Unchanged, except possibly to the extent of a small multiple of the
machine precision. (See
Further Comments.)
 2:
$\mathrm{x}\left({\mathbf{n}}\right)$ – double array

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.
 3:
$\mathrm{k}$ – int64int32nag_int scalar

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).
 4:
$\mathrm{el1n}$ – double scalar

The ${l}_{1}$norm (sum of absolute values) of the equation residuals.
 5:
$\mathrm{indx}\left({\mathbf{mpl}}\right)$ – int64int32nag_int array

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.
 6:
$\mathrm{ifail}$ – int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see
Error Indicators and Warnings).
Error Indicators and Warnings
Errors or warnings detected by the function:
 ${\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
nag_fit_glinc_l1sol (e02gb) again without changing the arguments.
 ${\mathbf{ifail}}=3$

The function has failed because of numerical difficulties; the problem is too illconditioned. Consider rescaling the unknowns.
 ${\mathbf{ifail}}=4$

Constraint: $\mathit{lde}\ge {\mathbf{n}}$.
Constraint: ${\mathbf{mpl}}\ge {\mathbf{m}}$.
Constraint: ${\mathbf{m}}\ge {\mathbf{n}}$.
Constraint: ${\mathbf{n}}\ge 2$.
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{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
Accuracy
The method is stable.
Further Comments
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 rankdeficient
$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.
Example
Suppose we wish to approximate in
$\left[0,1\right]$ a set of data by a curve of the form
which has nonnegative slope at the data points. Given points
$\left({t}_{i},{y}_{i}\right)$ we may form the equations
for
$\mathit{i}=1,2,\dots ,6$, for the
$6$ data points. The requirement of a nonnegative slope at the data points demands
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.)
Open in the MATLAB editor:
e02gb_example
function e02gb_example
fprintf('e02gb example results\n\n');
n = 4;
m = int64(6);
f = zeros(2*m,1);
e = zeros(n,2*m);
x = zeros(n,1);
f(1:m) = [0, 0.07, 0.07, 0.11, 0.27, 0.68];
for i = 1:m
t = 0.2*double(i1);
e(1:4,i) = [1,t,t^2,t^3];
e(1:4,m+i) = [0,1,2*t,3*t^2];
end
mxs = int64(100);
iprint = int64(0);
[e, x, k, l1nrm, indx, ifail] = e02gb(m, e, f, x, mxs, @monit, iprint);
function [] = monit(n, x, niter, k, elin)
fprintf('\n Results at iteration %d\n', niter);
fprintf('Constrained cubic polynomial:\n\n');
fprintf('p(t) = %7.4f + %7.4f t + %7.4f t^2 + %7.4f t^3\n\n',x(1:4));
fprintf('Norm of residuals = %12.5f\n', elin);
e02gb example results
Results at iteration 8
Constrained cubic polynomial:
p(t) = 0.0000 + 0.6943 t + 2.1482 t^2 + 2.1339 t^3
Norm of residuals = 0.00957
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015