e02gbf : NAG Library, Mark 26

e02gbf calculates an

l_{1}

solution to an over-determined system of linear equations, possibly subject to linear inequality constraints.

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 nagmk26.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)

Given a matrix

A

with

m

rows and

n

columns

(m \geq n)

and a vector

b

with

m

elements, the routine calculates an

l_{1}

solution to the over-determined system of equations

A x = 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) = \sum_{i = 1}^{m} |r_{i}|,

where the residuals

r_{i}

are given by

r_{i} = b_{i} - \sum_{j = 1}^{n} a_{i j} x_{j}, i = 1, 2, \dots, m .

Here

a_{i j}

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 (x)

subject to the set of inequality constraints

C x \geq 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

(t_{i}, y_{i})

is to be approximated by a linear combination of known functions

ϕ_{i} (t)

,

α_{1} ϕ_{1} (t) + α_{2} ϕ_{2} (t) + \dots + α_{n} ϕ_{n} (t),

in the

l_{1}

-norm, possibly subject to linear inequality constraints on the coefficients

α_{j}

of the form

C α \geq d

where

α

is the vector of the

α_{j}

and

C

and

d

are as in the previous paragraph. This is equivalent to finding an

l_{1}

solution to the over-determined system of equations

\sum_{j = 1}^{n} ϕ_{j} (t_{i}) α_{j} = y_{i}, i = 1, 2, \dots, m,

subject to

C α \geq d

.

Thus if, for each value of

i

and

j

, the element

a_{i j}

of the matrix

A

above is set equal to the value of

ϕ_{j} (t_{i})

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

α_{j}

. Note that the independent variable

t

above can, instead, be a vector of several independent variables (this includes the case where each of

ϕ_{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) - \sum_{i = 1}^{l} \min (0, c_{i}^{T} x - d_{i}),

where

γ

is a penalty parameter,

c_{i}^{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.

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

If on entry

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:

The method is stable.

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.

Suppose we wish to approximate in

[0, 1]

a set of data by a curve of the form

y = a x^{3} + b x^{2} + c x + d

which has non-negative slope at the data points. Given points

(t_{i}, y_{i})

we may form the equations

y_{i} = a t_{i}^{3} + b t_{i}^{2} + c t_{i} + d

for

i = 1, 2, \dots, 6

, for the

6

data points. The requirement of a non-negative slope at the data points demands

3 a t_{i}^{2} + 2 b t_{i} + c \geq 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.)

NAG Library Routine Document

e02gbf (glinc_l1sol)

▸▿ Contents

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results

NAG Library Routine Document

e02gbf (glinc_l1sol)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results