hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_fit_glinc_l1sol (e02gb)

Purpose

nag_fit_glinc_l1sol (e02gb) calculates an l1l1 solution to an over-determined 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 AA with mm rows and nn columns (mn)(mn) and a vector bb with mm elements, the function calculates an l1l1 solution to the over-determined system of equations
Ax = b.
Ax=b.
That is to say, it calculates a vector xx, with nn elements, which minimizes the l1l1-norm (the sum of the absolute values) of the residuals
m
r(x) = |ri|,
i = 1
r(x)=i=1m|ri|,
where the residuals riri are given by
n
ri = biaijxj,  i = 1,2,,m.
j = 1
ri=bi-j=1naijxj,  i=1,2,,m.
Here aijaij is the element in row ii and column jj of AA, bibi is the iith element of bb and xjxj the jjth element of xx.
If, in addition, a matrix CC with ll rows and nn columns and a vector dd with ll elements, are given, the vector xx computed by the function is such as to minimize the l1l1-norm r(x)r(x) subject to the set of inequality constraints CxdCxd.
The matrices AA and CC need not be of full rank.
Typically in applications to data fitting, data consisting of mm points with coordinates (ti,yi)(ti,yi) is to be approximated by a linear combination of known functions φi(t)ϕi(t),
α1φ1(t) + α2φ2(t) + + αnφn(t),
α1ϕ1(t)+α2ϕ2(t)++αnϕn(t),
in the l1l1-norm, possibly subject to linear inequality constraints on the coefficients αjαj of the form CαdCαd where αα is the vector of the αjαj and CC and dd are as in the previous paragraph. This is equivalent to finding an l1l1 solution to the over-determined system of equations
n
φj(ti)αj = yi,  i = 1,2,,m,
j = 1
j=1nϕj(ti)αj=yi,  i=1,2,,m,
subject to CαdCαd.
Thus if, for each value of ii and jj, the element aijaij of the matrix AA above is set equal to the value of φj(ti)ϕj(ti) and bibi is equal to yiyi and CC and dd are also supplied to the function, the solution vector xx will contain the required values of the αjαj. Note that the independent variable tt above can, instead, be a vector of several independent variables (this includes the case where each of φiϕ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
l
g(x) = γr(x)min (0, ciT xdi ),
i = 1
g(x) = γ r(x) - i=1 l min(0, ciT x-di ) ,
where γγ is a penalty parameter, ciT ciT  is the iith row of the matrix CC, and didi is the iith element of the vector dd. 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.

References

Bartels R H, Conn A R and Charalambous C (1976) Minimisation techniques for piecewise Differentiable functions – the ll 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 l1l1 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 l1l1 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

Parameters

Compulsory Input Parameters

1:     m – int64int32nag_int scalar
The number of equations in the over-determined system, mm (i.e., the number of rows of the matrix AA).
Constraint: m2m2.
2:     e(lde,mpl) – double array
lde, the first dimension of the array, must satisfy the constraint ldenlden.
The equation and constraint matrices stored in the following manner.
The first mm columns contain the mm rows of the matrix AA; element e(i,j)eij specifying the element ajiaji in the jjth row and iith column of AA (the coefficient of the iith unknown in the jjth equation), for i = 1,2,,ni=1,2,,n and j = 1,2,,mj=1,2,,m. The next ll columns contain the ll rows of the constraint matrix CC; element e(i,j + m)eij+m containing the element cjicji in the jjth row and iith column of CC (the coefficient of the iith unknown in the jjth constraint), for i = 1,2,,ni=1,2,,n and j = 1,2,,lj=1,2,,l.
3:     f(mpl) – double array
mpl, the dimension of the array, must satisfy the constraint mplmmplm.
f(i)fi, for i = 1,2,,mi=1,2,,m, must contain bibi (the iith element of the right-hand side vector of the over-determined system of equations) and f(m + i)fm+i, for i = 1,2,,li=1,2,,l, must contain didi (the iith element of the right-hand side vector of the constraints), where ll is the number of constraints.
4:     x(n) – double array
n, the dimension of the array, must satisfy the constraint n2n2.
x(i)xi must contain an estimate of the iith unknown, for i = 1,2,,ni=1,2,,n. If no better initial estimate for x(i)xi is available, set x(i) = 0.0xi=0.0.
5:     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 nn. If, on entry, mxs is zero or negative, the value returned by nag_machine_integer_max (x02bb) is used.
6:     monit – function handle or string containing name of m-file
monit can be used to print out the current values of any selection of its parameters. 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:     n – int64int32nag_int scalar
The number nn of unknowns (the number of columns of the matrix AA).
2:     x(n) – double array
The latest estimate of the unknowns.
3:     niter – int64int32nag_int scalar
The number of iterations so far carried out.
4:     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:     el1n – double scalar
The l1l1-norm of the current residuals of the over-determined system of equations.
7:     iprint – int64int32nag_int scalar
The frequency of iteration print out.
iprint > 0iprint>0
monit is called every iprint iterations and at the solution.
iprint = 0iprint=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:     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, nn (the number of columns of the matrix AA).
Constraint: n2n2.
2:     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 + lm+l, where ll is the number of constraints (which may be zero).
Constraint: mplmmplm.

Input Parameters Omitted from the MATLAB Interface

lde w iw

Output Parameters

1:     e(lde,mpl) – double array
ldenlden.
Unchanged, except possibly to the extent of a small multiple of the machine precision. (See Section [Further Comments].)
2:     x(n) – double array
The latest estimate of the iith unknown, for i = 1,2,,ni=1,2,,n. If ifail = 0ifail=0 on exit, these are the solution values.
3:     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:     el1n – double scalar
The l1l1-norm (sum of absolute values) of the equation residuals.
5:     indx(mpl) – int64int32nag_int array
Specifies which columns of e relate to the inactive equations and constraints. indx(1)indx1 up to indx(k)indxk number the active columns and indx(k + 1)indxk+1 up to indx(mpl)indxmpl number the inactive columns.
6:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
The constraints cannot all be satisfied simultaneously: they are not compatible with one another. Hence no solution is possible.
  ifail = 2ifail=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 parameters.
  ifail = 3ifail=3
The function has failed because of numerical difficulties; the problem is too ill-conditioned. Consider rescaling the unknowns.
  ifail = 4ifail=4
On entry, one or more of the following conditions are violated:
  • mn2mn2,
  • or mplmmplm,
  • or iw3 × mpl + 5 × n + n2 + (n + 1) × (n + 2) / 2iw3×mpl+5×n+n2+(n+1)×(n+2)/2,
  • or ldenlden.
Alternatively elements 11 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 AA or CC.

Accuracy

The method is stable.

Further Comments

The effect of mm and nn 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 nn and the total time taken is approximately proportional to mn2mn2.
Linear dependencies among the rows or columns of AA and CC are not necessarily a problem to the algorithm. Solutions can be obtained from rank-deficient AA and CC. 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

function nag_fit_glinc_l1sol_example
m = int64(6);
e = [1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0;
     0, 0.2, 0.4, 0.6, 0.8, 1, 1, 1, 1, 1, 1, 1;
     0, 0.04, 0.16, 0.36, 0.64, 1, 0, 0.4, 0.8, 1.2, 1.6, 2;
     0, 0.008, 0.064, 0.216, 0.512, 1, 0, 0.12, 0.48, 1.08, 1.92, 3];
f = [0;
     0.07;
     0.07;
     0.11;
     0.27;
     0.68;
     0;
     0;
     0;
     0;
     0;
     0];
x = [0;
     0;
     0;
     0];
mxs = int64(50);
iprint = int64(0);
[eOut, xOut, k, el1n, indx, ifail] = nag_fit_glinc_l1sol(m, e, f, x, mxs, @monit, iprint)

function [] = monit(n, x, niter, k, elin)

  fprintf('\n Results at iteration %d\n', niter);
  fprintf('X-Values\n');
  disp(transpose(x));
  fprintf('Norm of residuals = %12.5f\n', elin);
 

 Results at iteration 8
X-Values
         0    0.6943   -2.1482    2.1339

Norm of residuals =      0.00957

eOut =

  Columns 1 through 9

    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000         0         0         0
         0    0.2000    0.4000    0.6000    0.8000    1.0000    1.0000    1.0000    1.0000
         0    0.0400    0.1600    0.3600    0.6400    1.0000         0    0.4000    0.8000
         0    0.0080    0.0640    0.2160    0.5120    1.0000         0    0.1200    0.4800

  Columns 10 through 12

    0.0000   -0.0000    0.0000
    1.0000    1.0000    1.0000
    1.2000    1.6000    2.0000
    1.0800    1.9200    3.0000


xOut =

         0
    0.6943
   -2.1482
    2.1339


k =

                    4


el1n =

    0.0096


indx =

                    6
                    2
                    9
                    1
                    5
                   11
                   10
                    4
                    3
                    7
                   12
                    8


ifail =

                    0


function e02gb_example
m = int64(6);
e = [1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0;
     0, 0.2, 0.4, 0.6, 0.8, 1, 1, 1, 1, 1, 1, 1;
     0, 0.04, 0.16, 0.36, 0.64, 1, 0, 0.4, 0.8, 1.2, 1.6, 2;
     0, 0.008, 0.064, 0.216, 0.512, 1, 0, 0.12, 0.48, 1.08, 1.92, 3];
f = [0;
     0.07;
     0.07;
     0.11;
     0.27;
     0.68;
     0;
     0;
     0;
     0;
     0;
     0];
x = [0;
     0;
     0;
     0];
mxs = int64(50);
iprint = int64(0);
[eOut, xOut, k, el1n, 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('X-Values\n');
  disp(transpose(x));
  fprintf('Norm of residuals = %12.5f\n', elin);
 

 Results at iteration 8
X-Values
         0    0.6943   -2.1482    2.1339

Norm of residuals =      0.00957

eOut =

  Columns 1 through 9

    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000         0         0         0
         0    0.2000    0.4000    0.6000    0.8000    1.0000    1.0000    1.0000    1.0000
         0    0.0400    0.1600    0.3600    0.6400    1.0000         0    0.4000    0.8000
         0    0.0080    0.0640    0.2160    0.5120    1.0000         0    0.1200    0.4800

  Columns 10 through 12

    0.0000   -0.0000    0.0000
    1.0000    1.0000    1.0000
    1.2000    1.6000    2.0000
    1.0800    1.9200    3.0000


xOut =

         0
    0.6943
   -2.1482
    2.1339


k =

                    4


el1n =

    0.0096


indx =

                    6
                    2
                    9
                    1
                    5
                   11
                   10
                    4
                    3
                    7
                   12
                    8


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013